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18  Subject  terms  (Continue  on  reverse  if  necessary  and  identify  by  bloch  number) 

Regional  wave  propagation  Scattering  P  coda 

Body  waves  Norres 

Surface  waves  Finite  difference 


')  ABS'RACT  iConfinue  on  ceverse  .E  necessary  jnd  idfAfify  by  6(oc*  number)  Broad— band,  digital,  tcleselsmlc  P— WaVe 
data  from  the  University  of  California,  Berkeley,  stations  BKS,  MHC,  and  SAO  are  used 
to  study  the  character  of  teleselsmlc  P  coda.  Velocity  heterogeneity  in  structure 
near  these  stations  gives  rise  to  a  variety  of  scattering  effects  which  include 
azimuthly  dependent  receiver  functions,  large  tangential  component  particle  motions, 
and  high  levels  of  P  coda.  P  codas  for  shallow  and  deep  teleseisms  are  analyzed  for 
three  frequency  bands  centered  at  1.5,  0.75,  and  0.375  hz.  Observed  coda  decay  rate 
and  level  are  fit  with  an  energy  flux  theory  to  parameterize  scattering  effects  near  the 
source  and  at  the  receiver.  Coda  decay  rate  and  level  for  deep  earthquakes  are  seen  to 
be  stable  observables  between  stations,  implying  that  the  scattered  coda  field  is  homo¬ 
geneously  distributed  over  the  receiver  network  validating  a  major  assumption  In  the 
energy  flux  model.  Scattering  Q,  is  not  particularly  well  constrained  by  the  data 
and  is  inferred  to  be  affected  by  focusing  and  defocuslng  effects  of  the  direct  P-waves. 
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There  is  a  pronounced  difference  in  coda  decay  rate  between  shallow  and  deep  earth¬ 
quakes  for  1.5  hz  coda  where  coda  decay  rate  for  deep  events  is  about  a  factor  of 
two  greater  than  for  shallow  events.  There  is  no  difference  in  coda  decay  rate  be¬ 
tween  the  two  populations  of  events  in  the  lower  frequency  bands.  The  observation 
at  1.5  hz  is  consistent  with  a  simple  model  of  P-wave  coda  formation  where  most  coda 
forms  in  a  scattering  layer  near  the  surface  of  the  earth  in  regions  near  the  source 
and  near  the  receiver.  The  homogeneity  of  coda  decay  at  0.75  and  0.375  hz  for  deep 
shallow  events  suggests  that  a  large,  near-receiver  scattering  mechanism,  such  as 
surface  wave  formation,  dominates  the  coda. 

Events  recorded  by  NORESS  are  analyzed  to  determine  the  velocities  and  times  of 
coherent  crustal  phases  to  be  used  in  accurate  locations  and  determination  of  crustal 
structure.  The  events  straddle  two  azimuthal  directions,  and  range  in  distance  from 
50  to  300  km,  thus  sampling  two  cross-sections  of  the  crust  and  upper  mantle;  one  in 
the  Caledonides  to  the  north  west  of  NORESS,  and  the  other  to  the  south  along  the 
eastern  rim  of  the  Oslo  Graben.  F-k  analysis  is  used  to  search  for  coherent  phases 
in  the  data.  All  vertical  channels  are  then  stacked  with  the  appropriate  time 
delays  calculated  from  the  phase  velocities  and  azimuths  obtained  for  the  coherent 
arrivals.  Finally,  a  composite  seismogram  is  made  for  each  event  by  piecing  together 
time  sections  from  the  stacks,  where  each  time  piece  represents  a  coherent  arrival. 
All  detectable  coherent  phases  are  enhanced  in  the  composite  seismograms  and  on 
closer  inspection  it  is  possible  to  improve  some  of  the  locations  from  the  NORESS 
bulletin.  This  causes  a  number  of  events  to  group  together  to  the  same  location. 
Travel  time  curves  are  constructed  from  plane  layered  homogeneous  models  to  match 
travel  times  and  phase  velocities  in  the  composite  seismograms. 

Numerical  simulations  of  teleselsmlc  wave  propagation  in  a  heterogeneous  layer  over 
a  homogeneous  half-space  are  conducted  to  further  understanding  of  teleselsmlc  coda 
generation  and  complex  teleselsmlc  particle  motions.  Acoustic  and  elastic  finite 
difference  synthetics  are  generated  for  more  than  150  different  layer  over  a  half¬ 
space  models.  The  models  vary  in  scattering  layer  thickness  (L) ,  heterogeneity 
correlation  distance(s)  (a,  or  a  ,  a  ) ,  and  heterogeneity  standard  deviation  (o). 

The  synthetic  data  are  examined  using  relative  coda  intensity  envelopes,  to  determ¬ 
ine  relationships  that  exist  between  the  model  parameters  and  the  level  and  rate  of 
coda  decay,  and  frequency-wavenumber  analysis  to  determine  the  wave  types  present  in 
the  scattered  wavefield.  The  level  of  scattered  energy  is  found  to  be  controlled 
by  the  ka  value  and  the  magnitude  of  the  standard  deviation  of  the  perturbations. 

Coda  levels  increase  with  ka  from  ka<l  to  ka~  1,  and  decrease  with  increasing 
ka  for  ka  >1.  Scattered  energy  levels  are  found  to  always  increase  with  increasing 
A  more  complex  relationship  is  found  to  exist  between  the  coda  level  and  decay  rate, 
and  variations  in  the  scattering  layer  thickness.  For  models  that  vary  in  scatter¬ 
ing  layer  thickness  alone,  coda  level  and  decay  rate  changes  did  not  occur  independ¬ 
ently  of  the  values  of  ka  and/or  O  .  The  rate  of  coda  decay  is  found  to  be  controll¬ 
ed  by  the  aspect  ratio  (a  /a  )  of  the  heterogeneities  (spatial  anisotropy).  The  two 
extreme  models  are  those  wltfi  homogeneous  plane-layers,  %rtiich  exhibit  the  most  rapid 
coda  decay  rates,  and  those  with  spatially  isotropic  heterogeneities  (a  /a  =1), 
which  exhibit  the  slowest  decay  rates.  Acoustic  and  elastic  models  exhibif  similar 
coda  intensity  envelope  characteristics.  Apparent  scattering  attenuation  of  the 
direct  pulse  is  found  to  be  strongest  for  models  with  spatially  isotropic  heterogen¬ 
eities.  Frequency-wavenumber  analysis  of  the  elastic  synthetics  reveals  the  mechan¬ 
ism  responsible  for  the  generation  of  coda  in  the  synthetic  data.  For  models  with 
spatially  isotropic  heterogeneities,  the  coda  is  composed  largely  of  low  apparent 
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velocity  (LAV)  waves  (scattered  S  and  possibly  body-to-surf ace  scattered  waves) 
%rtjile  the  coda  for  models  with  homogeneous  plane-layers  is  composed  solely  of 
vertically  propagating  layer  reverberations.  Decreasing  the  heterogeneity  aspect 
ratio  increases  the  total  amount  of  scattered  LAV  energy  produced.  The  onset  time 
of  the  LAV  energy  is  also  controlled  by  the  heterogeneity  aspect  ratio.  For  models 
with  spatially  Isotropic  heterogeneities,  the  LAV  energy  appears  immediately  after 
the  initial  arrival.  For  models  with  homogeneous,  plane-layers,  LAV  energy 
never  appears. 

Teleseismlc  data  recorded  at  the  NORESS  seismic  array  are  examined  to  determine  the 
degree  and  scale  of  lithospheric  heterogeneities  beneath  the  array,  and  to  locate 
local  scattering  sources.  The  data  are  examined  using  relative  coda  Intensity 
envelopes,  frequency-wavenumber  spectra,  and  particle  motion  polarizations.  The 
coda  intensity  envelopes  do  not  exhibit  any  clear  ka  trends;  the  0.5  Hz  pass-band 
data  does,  however,  contain  higher  energy  levels  .and  a  high  energy  packet  ~  15 
seconds  after  the  first  arrival.  Frequency-wavenumber  analysis  shows  the  most  ot 
the  teleseismlc  coda  is  vertically  propagating.  The  window  centered  near  15 
seconds,  however,  contains  energy  travelling  to  the  north-east  at  a  low  apparent 
velocity.  Particle  motions  in  this  window  follow  a  retrograde  elliptical  pattern 
characteristic  of  Rayleigh  surface  waves.  The  scattering  source  presumed  to  be 
responsible  for  the  observed  body-to-surface  wave  scattered  energy  is  Lake  Mjosa, 
located  approximately  50  km  west-southwest  of  the  array. 
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Report  Summary 


Task  Objectives 


Characterization  of  local  and  regional  wave  propagation  in  the  lithosphere 
around  the  NORESS  array. 

Theoretical  and  observational  study  of  teleseismic  P-wave  coda  near  seismic 
receivers. 

Technical  Problem 

Small  aperture  seismic  arrays,  like  the  NORESS  array,  are  central  in  the 
problem  of  underground  nuclear  test  detection,  discrimination  and  yield 
estimation.  Array  geometries  allow  an  analysis  of  the  wave  field  of 
propagating  high  frequency  regional  waves  which  can  be  used  to  determine 
source  location  and  size.  Integral  to  the  location  problem  is  detailed 
knowledge  of  wave  propagation  through  the  crust  and  mantle  between  the 
source  and  array.  The  accuracy  of  a  location  is  intimately  tied  to  simple 
identification  of  discrete  phases  which  arrive  at  the  array  which  are  used  to 
estimate  source  distance  and  azimuth.  This  work  addressed  several 
fundamental  seismological  problems  in  wave  propagation  dealing  with  the 
identification  of  important  phases  in  the  NORESS  data  and  propagation  of 
those  phases  through  the  crust  and  mantle. 

We  examined  two  broad  areas  of  wave  propagation  important  in 
characterizing  the  nature  of  wave  propagation  at  regional  distances.  The  first 
was  a  fundamental  observational  study  of  seismic  wave  propagation  near 
NORESS.  Basic  observational  data  were  needed  to  analyze  the  excitation  and 
behavior  of  common  regional  phases  such  as  Pn,  Pg,  Sn,  Lg  and  Rg.  The 
other  general  area  approached  the  problem  of  the  effect  of  wave  scattering  in 
the  crust  and  mantle  and  how  it  could  be  characterized  using  a  "simpler"  data 
set  consisting  of  teleseismic  P  waves.  A  working  assumption  was  that 
observation  of  near-receiver  scattering  in  teleseismic  P  wave  coda  could  help 
characterize  the  general  nature  of  scattering  due  to  near-receiver  velocity 


heterogeneity  for  local  and  regional  events.  Teleseismic  NORESS  data  were 
incorporated  in  this  aspect  of  the  observational  study. 

General  Methodology 

The  study  of  local  and  regional  phases  was  performed  through  direct 
observation  using  NORESS  data  and  applying  array  processing  techniques. 
Theoretical  studies  of  wave  propagation  in  layered  media  were  performed 
using  computer  simulation  to  examine  the  expected  nature  of  wave 
propagation  in  simple  crustal  models.  Finite  Difference  simulations  were 
also  performed  to  examine  the  problem  of  scattering  in  a  layered  halfspace 
assuming  incident  teleseismic  P  plane  waves.  Teleseismic  data  at  NORESS 
were  examined  to  infer  the  nature  of  scattering  in  the  lithosphere  under 
NORESS. 

Technical  Results 

Section  1  describes  the  results  of  careful  examination  of  regional  phases  in  the 
NORESS  data  set.  Array  techniques  were  used  to  decompose  the  seismogram 
into  constituent  phases  for  analysis  in  the  time  domain.  Gross  designation  of 
regional  phases  are  commonly  used  in  routine  seismic  analysis.  Phases  such 
as  Pg  and  Lg  are  picked  from  the  data  using  rough  estimates  of  phase  velocity 
and  arrival  time  and  are  then  incorporated  into  location  algorithms.  Our 
study  showed  that  the  crust  in  Fennoscandia  generates  a  suite  of  clear, 
understandable  arrivals  which  are  poorly  described  by  previous  simple  labels. 
These  P,  S,  and  Rayleigh  wave  arrivals  contain  important  information  on 
phase  velocity  and  timing  which  can  be  used  for  more  accurate  event  distance 
estimates.  Indeed,  we  found  that  standard  procedures  used  to  generate  the 
NORESS  catalog  often  misidentified  S  and  P  phases  resulting  in  mislocations 
of  80  km  or  more.  We  also  found  that  much  of  the  wave  propagation  could 
be  understood  with  fairly  simple  layered  models.  The  effect  of  scattering  due 
to  heterogeneity,  while  important  for  some  azimuths  and  ranges,  was  not 
overwhelming.  Thus,  future  array  observations  n\ay  be  able  to  lake 
advantage  of  more  accurate  processing  to  identify  phases  in  the  regional  wave 
form. 


Sections  2,  3,  and  4  contain  results  on  the  nature  of  velocity  heterogeneity 
seen  in  the  lithosphere  near  receivers  from  incident  teleseismic  P  waves. 
This  work  grew  out  of  studies  performed  under  the  previous  contract  on  the 
nature  of  scattering  which  produces  teleseismic  P  wave  coda.  In  Section  2,  an 
observational  study  was  performed  to  examine  two  assumptions  in  the 
energy  flux  model  of  coda  generation  developed  in  past  work.  Data  from  the 
Berkeley  broad-band  network  were  analyzed  to  determine  if  coda  levels  and 
decays  are  constant  over  a  large  region.  This  would  substantiate  the 
assumption  of  homogeneity  of  coda  energy  distribution  behind  the  incident 
wavefront,  important  in  developing  the  energy  flux  theory.  We  found  that 
this  was  indeed  true.  However,  we  found  that  assuming  that  the  direct  wave 
is  little  affected  by  heterogeneity  is  not  very  good.  Near-receiver  effects  serve 
to  focus  and  defocus  the  direct  wave  invalidating  the  normalization  required 
by  the  theory  to  infer  the  scattering  Q  parameter.  We  also  found  that 
teleseismic  P  wave  coda  can  contain  a  significant  portion  due  to  scattering 
near  the  source. 

Section  3  describes  an  extensive  theoretical  study  of  teleseismic  P  wave  coda 
using  2D  finite  difference  techniques  to  solve  the  acoustic  and  elastic  wave 
equation.  A  scattering  layer  over  an  isotropic,  non-scattering  halfspace  was 
assumed.  We  were  interested  in  determining  which  factors  were  most 
important  in  controlling  coda  decay  and  level.  Langston(1989)  asserted  that 
the  distribution  of  scatterers  in  a  layer  could  be  very  important  for  both  the 
level  and  decay  of  coda.  The  simulations  showed  a  surprising  result  that 
isotropic  scattering  produced  the  same  coda  decay  for  wide  variations  in  layer 
thickness  and  amplitude  of  the  velocity  perturbations.  The  rate  of  coda  decay 
was  most  sensitive  to  the  anisotropy  of  the  scatterers.  Models  which  were 
plane-layered-like  produced  the  highest  decay  rates.  Decay  rates  were  slowest 
for  the  isotropic  scatterers.  This  presents  the  important  possibility  that 
reconnaissance  of  the  gross  geometry  of  scatterers  in  the  lithosphere  may  be 
possible  by  examining  coda  decay  rates  and  the  directionality  of  the  scattered 
coda  field. 

This  was  attempted  using  teleseismic  P  waves  recorded  by  NORHSS  and  is 
described  in  Section  4.  It  was  found  that  teleseismic  coda  waves  from  deep 
events  are  generally  dominated  by  propagation  in  near-vertical  directions, 
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implying  that  local  scattering  is  probably  dominated  by  plane-layered-like 
structure.  However,  a  discrete  and  significant  P  to  Rg  scattering  event  was 
found  in  the  data  and  was  associated  with  the  Lake  Mjosa  topographic  feature 
near  the  NORESS  array. 

Important  Findings  and  Conclusions 

Wave  propagation  in  Fennoscandia  is  generally  understandable  in  terms  of 
simple  plane  layered  structure  models.  Although  the  data  show  deviations 
from  such  models  in  terms  of  effects  of  layer  dip,  discrete  phases  are  generated 
by  these  crust  and  mantle  structures  and  are  not  significantly  degraded  by 
heterogeneity.  Phase  attributes,  such  as  time  and  phase  velocity,  can  be  used 
to  accurately  identify  arrivals.  It  was  found  that  S  waves  are  often 
misidentified  at  NORESS  resulting  in  poor  event  locations. 

The  nature  of  lithospheric  heterogeneity  which  produces  P  wave  coda  may  be 
amenable  to  analysis  since  it  was  found  that  scatterer  anisotropy  appeared  to 
be  the  strongest  effect  in  controlling  coda  decay  rate  in  the  wave  propagation 
simulations.  The  NORESS  data  suggest  that  scattering  in  the  vicinity  of  the 
array  is  caused  by  highly  anisotropic  scatterers  (plane-layered-like)  with  the 
exception  of  a  discrete  scattering  events  caused  by  P-to-Rg  scattering  off  of 
topographic  features. 

Significant  Hardware  Development 

N/A 

Special  Comments 
N/A 

Implications  for  Further  Research 

The  careful  analysis  of  a  commonly  available  NORESS  data  set  points  out  that 
there  is  still  much  to  learned  about  the  propagation  of  regional  phases 
through  direct  observation.  Other  arrays  of  the  area  are  currently  being 


studied  in  a  companion  project  to  broaden  our  experience  with  the 
generation  of  high-frequency  regional  phases.  A  challenge  will  be  to 
implement  these  results  into  fast,  array  algorithms  to  refined  location 
estimates. 
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Abstract 

Events  recorded  by  P^ORESS  are  analyzed  to  determine  the  velocities  and 
times  of  coherent  crustal  phases  to  be  used  in  accurate  locations  and  determination 
of  crustal  structure.  The  events  straddle  two  azimuthal  directions,  and  range  in  dis¬ 
tance  from  50  to  300  km,  thus  sampling  two  cross-sections  of  the  crust  and  upper 
mantle;  one  in  the  Caledonides  to  the  north  west  of  NORESS,  and  the  other  to  the 
south  along  the  eastern  rim  of  the  Oslo  Graben.  F-k  analysis  is  used  to  search  for 
coherent  phases  in  the  data.  All  vertical  channels  are  slant-stacked  using  the  phase 
velocities  and  azimuths  obtained  in  the  windows  containing  coherent  arrivals.  A 
composite  seismogram  is  then  made  for  each  event  by  piecing  together  slant-stacked 
sections  representing  all  the  detectable  pha.ses.  Coherent  phases  are  greatly 
enhanced  by  the  procedure,  and  on  closer  inspection  it  is  possible  to  improve  some 
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of  the  locations  from  the  NORESS  bulletin.  This  causes  a  number  of  events  to 
group  together  to  the  same  location.  PmP  and  SmS  are  the  largest  arrivals  on  all 
records  from  distances  greater  than  80  km.  Closer  in,  the  largest  arrivals  are  P  and 
S  waves  turning  or  reflecting  in  the  upper  crust  Rg  waves  are  also  prominent  in 
records  from  that  distance  range.  Pn  is  detected  at  distances  around  200  km  and 
further  and  has  a  very  high  phase  velocity  in  the  Caledonian  record  section, 
between  9  and  10  km/s.  Sn  can  also  be  found  in  some  events.  Double  crustal 
reverberations  of  P  and  S  are  observed  at  distances  greater  than  200  km.  Travel 
time  curves  are  constructed  from  plane  layered  homogeneous  models  to  match 
travel  times  and  phase  velocities  in  the  composite  seismograms. 

Introduction 

Studies  of  wave  propagation  from  regional  seismic  events  and  seismic  refrac¬ 
tion  lines  have  generally  been  hampered  by  low  signal-to-noise  ratio  and  by  poor 
phase  identification.  With  the  use  of  high-frequency,  small-aperture  seismic  arrays, 
such  studies  can  be  greatly  improved.  Array  analysis  not  only  enhances  the  s/n 
ratio  to  bring  out  more  details  in  the  seismograms,  but  also  allows  the  association 
of  phase  velocities  and  backazimuths  with  the  arrivals.  At  regional  distances,  the 
prominent  phases  in  the  seismograms  are  the  Pg  and  Lg  wave  trains,  which  are 
made  up  of  a  sequence  of  multiply  reflected  postcritical  P  and  S  waves  trapped  in 
the  crustal  waveguide.  The  main  reflectors  are  the  free  surface  and  the  Moho,  but 
internal  discontinuities  in  the  crust  also  act  to  increase  the  complexity  of  the  wave 
trains.  The  number  of  multiples  increases  with  source  distance  as  the  critical  dis¬ 
tance  for  each  multiple  reflection  is  reached  (Bouchon,  1982).  The  relative  arrival 
times  and  amplitudes  of  the  phases  in  the  wave  trains  are  sensitive  to  crustal 


structure  and  source  depth.  They  are  therefore  ideally  suited  for  studying  crustal 
structure  and  source  effects  in  regional  events. 

The  present  study  uses  regional  events  recorded  by  the  NORESS  array,  which 
was  designed  to  enhance  the  detection  of  such  events  (Mykkeltveit  et  al.,  1983),  to 
study  structure  and  its  effects  on  wave  propagation  at  regional  distances.  The  array 
is  located  in  the  Baltic  Shield,  northeast  of  the  Oslo  Graben,  which  was  fonned  by 
a  continental  rifting  episode  during  the  Permian.  Volcanics  and  intrusions  are  pre¬ 
valent  in  the  Graben  area,  and  results  from  previous  seismic  and  gravity  studies 
indicate  undulations  of  several  kilometers  in  the  Moho  under  the  Oslo  Graben 
(Ramberg,  1976;  Wessel  and  Husebye,  1987;  Kanestrom  and  Haugland,  1971). 

The  events  were  chosen  from  the  NORESS  bulletin  based  on  location  and 
local  magnitude  to  allow  sampling  of  the  crust  on  two  different  profiles  extending 
from  the  array.  The  location  of  events,  which  range  in  magnitude  from  1.7  to  3.8, 
is  shown  in  Figure  I.  The  NW-SE  crustal  profile,  delineated  by  the  event  distribu¬ 
tion  extends  into  the  Caledonian  mountain  range,  while  the  N-S  profile,  runs  along 
the  eastern  rim  of  the  Oslo  Graben.  The  events  are  subject  to  f-k  analysis  to  search 
for  coherent  arrivals;  these  are  then  used  to  infer  crustal  structure  and  to  improve 
event  locations,  where  necessary.  The  seismograms  are  compared  with  travel  time 
curves  calculated  for  two  different  crustal  models  obtained  from  previous  refraction 
studies  in  the  region  (Sellevoll  and  Warrick,  1971;  Kanestrom  and  Haugland,  1971). 

Data  Analysis 

In  order  to  find  coherent  phases  in  the  data  and  to  identify  their  phase  veloci¬ 
ties.  f-k  analysis  is  performed  in  sliding  time  windows  on  the  25  vertical  com¬ 
ponents  of  the  array.  Window  lengths  and  their  locations  in  time  are  chosen  after 
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inspection  of  the  record  sections,  so  that  the  window  will  contain  the  entire  phase, 
which  is  at  the  center  of  the  window.  The  data  are  then  multiplied  by  a  Hanning 
smoothing  function  and  Fourier  transformed  in  time  to  find  the  spectral  peaks 
within  the  window.  Two-dimensional  Fourier  transform  on  the  space  dimensions  is 
then  calculated  at  the  peak  frequencies  to  obtain  the  power  spectrum  in  the 
wavenumber  plane.  If  the  window  contains  a  coherent  arrival,  then  there  will  be  a 
distinct  maximum  peak  of  power  somewhere  in  the  wavenumber  plane.  The  loca¬ 
tion  of  the  peak  determines  the  phase  velocity  and  the  direction  from  which  the 
wave  is  coming.  Phase  velocities  and  backazimuths  obtained  for  the  various  fre¬ 
quency  peaks  within  a  window  often  show  some  variation.  This  can  be  due  to 
interference  between  phases  with  different  travel  paths,  caused  by  layering  or 
heterogeneities  in  the  crust.  Lack  of  resolution  can  also  be  the  cause  if  phase  velo¬ 
cities  are  high.  In  this  situation  the  phase  velocity  and  backazimuth  of  the  fre¬ 
quency  with  the  greatest  px)wer  is  chosen  to  represent  the  arrival.  When  there  is 
variation  in  backazimuth  between  the  main  P  and  S  phases  within  one  event,  the 
backazimuth  of  the  largest  phase  is  chosen  to  represent  the  event. 

With  the  backazimuth  and  phase  velocities  for  all  detectable  phases  of  an 
event  in  hand,  the  vertical  channels  are  slant  stacked  and  corrected  to  zero  level. 
Time  sections  from  the  different  slant  stacks  are  then  pieced  together  to  form  one 
composite  seismogram  as  shown  in  Figure  2.  The  composite  has  a  greatly  improved 
signal-to-noise  ratio  so  that  the  major  crustal  phases  clearly  stand  out  and  some  of 
the  more  obscure  ones,  such  as  Pn,  become  apparent.  In  some  of  the  more  distant 
events  with  small  magnitudes,  the  signal-to-noise  ratio  is  only  2  -  5.  In  such  cases, 
the  f-k  analysis  is  performed  on  the  raw  data,  but  the  channels  are  high-pass  filtered 
with  a  3-pole,  zero-phase  Butterworth  filter  before  slant  stacking  to  eliminate  the 


low-frequency  noise. 

When  the  channels  are  slant-stacked,  coherent  arrivals  with  the  preferred 
backazimuth  and  incidence  angle  are  enhanced,  while  the  incoherent  amplitude,  and 
coherent  amplitude  arriving  with  a  different  backarimuth  and  incidence  angle,  are 
diminished.  If  a  substantial  variation  in  backazimuth  and  phase  velocity  occurs 
within  an  arrival  the  event  will  appear  to  have  low  coherency.  For  a  quantitative 
estimate  of  the  coherency  in  each  event,  and  in  order  to  compare  the  effects  of 
different  travel  paths,  the  envelopes  of  ail  vertical  channels  are  slant-stacked  and 
time  sections  pieced  together  in  the  same  manner  as  the  seismograms;  this  compo¬ 
site  envelope  approximates  the  true  total  amplitude.  The  envelope  of  the  composite 
seismogram,  which  approximates  the  true  coherent  amplitude  is  then  subtracted  to 
obtain  the  incoherent  amplitude. 


Event  Locations 

After  slant  stacking  the  events  presented  in  Figure  1,  it  became  apparent  that 
several  were  mislocated,  some  by  as  much  as  80  km.  The  composite  seismograms 
for  the  events  in  question  are  presented  in  Figures  3,  4  and  5,  where  the  arrows 
indicate  the  time  picks  of  coherent  arrivals  from  the  NORESS  bulletin.  Some  of 
the  events  were  recorded  during  a  period  of  unknown  time  drift  of  <  2  sec,  occur¬ 
ring  in  the  on-line  data.  This  time  drift  is  indicated  in  the  figures  by  the  time  bars. 
The  two  thicker  arrows  in  each  event  indicate  the  arrivals  of  the  P  and  Lg  waves, 
which  were  used  to  locate  the  event. 

There  are  two  main  reasons  for  the  mislocations.  The  first  main  reason 
involves  the  picking  of  the  largest  S  wave  arrival  on  the  record  as  the  arrival  of  the 
Lg  wave.  The  two  events  in  Figure  2  have  two  main  S  arrivals:  Sg,  which  is  used 
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here  to  represent  an  S  wave  that  turns  or  reflects  in  the  crust,  with  a  phase  velocity 
of  3.7  km/s,  and  the  longer  period  Rg  surface  wave,  with  a  velocity  of  3  km/s.  A 
schematic  plot  explaining  the  travel  paths  of  the  various  phases  is  shown  in  Figure 
6.  By  picking  the  larger  amplitude  Rg  wave  as  the  onset  of  Lg  in  event  88316,  the 
source  distance  is  overestimated  by  40  km.  Figure  3  shows  the  same  misinterpreta¬ 
tion  of  the  onset  of  Lg  for  four  nearly  identical  events,  with  three  major  S  wave 
arrivals:  Sg,  the  Moho  reflected  SmS,  and  Rg.  Depending  on  which  phase  is  picked 
as  the  Lg  arrival,  the  events  are  misplaced  from  the  actual  distance  of  about  85  km, 
by  as  much  as  80  km. 

The  second  main  reason  for  mislocation  is  demonstrated  in  Figure  5,  which 
shows  four  events  with  the  same  relative  arrival  times  but  located  over  a  70  km  dis¬ 
tance  range.  The  two  closer  events,  882205  and  88286,  have  local  magnitudes  of 
2.3  and  2.2,  whereas  the  other  two,  88257  and  87294  have  local  magnitudes  2.6  and 
2.8.  Lg  is  picked  at  about  the  same  time  in  all  events,  at  the  arrival  of  SmS.  Pn  is 
picked  as  the  first  P  arrival  in  the  two  more  distant  events,  but  due  to  its  smaller 
amplitude  it  is  missed  in  events  88286  and  88205;  2xFhnP  is  picked  for  event  88205 
and  the  time  pick  for  88286  is  late  in  the  PmP  phase.  Thus  small  distant  events 
tend  to  be  mislocated  closer  than  they  actually  are  when  the  first  arriving  P  phases 
are  not  detected.  The  actual  distance  for  the  four  events  is  about  240  km,  depending 
on  velocity  model.  The  mislocated  events  show  a  distribution  in  backazimuth,  of  up 
to  6  degrees,  indicating  some  errors  in  backazimuth  also. 

Figure  1  shows  three  events  with  local  magnitudes  above  3.  These  are  errone¬ 
ous  magnitudes,  and  when  recalculated  according  to  Bath  (1976)  they  decreased  to 
between  1.1  and  1.8.  NORESS  magnitudes  were  used  for  all  the  other  events.  The 
events  are  replotted  in  Figure  7  in  their  relocated  positions. 
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Record  Sections 


NW-SE  Section 

The  NW-SE  record  section  has  five  events,  which  range  in  distance  and  backa- 
zimuth  from  about  50  to  300  km,  and  293  to  314  degrees,  respectively.  Two 
events,  87183  and  85313,  have  been  high-pass  filtered  for  clearer  phase 
identification.  The  record  section  is  shown  in  Figures  8  and  9.  Pg  is  the  first 
arrival  in  the  two  closest  events,  but  PmP  takes  over  before  the  source  distance  of 
the  third  event  is  reached.  Event  88244,  which  is  the  closest  at  about  50  km  dis¬ 
tance,  has  well-defined  Pg  and  Sg  phases  and  a  large  Rg  phase,  but  no  mantle 
reflection  is  apparent.  Event  88133,  on  the  other  hand,  is  close  to  the  critical  dis¬ 
tance  for  mantle  reflection  and  therefore  has  large  PmP  and  SmS  arrivals  following 
Pg  and  Sg.  The  Rg  phase  is  also  still  prominent. 

Pn  is  the  first  arrival  on  the  remaining  events  with  a  phase  velocity  of  about 
10  km/s.  Sn  is  also  found  in  the  two  farthest  events  with  a  phase  velocity  of  about 
5.2  km/s.  A  large  PmP  phase  follows  Pn  in' all  three  events,  with  a  phase  velocity 
of  about  7.2  km/s,  and  2xPmP  with  a  phase  velocity  of  around  8  km/s  starts  build¬ 
ing  up  from  a  distance  of  200  km  and  further.  The  corresponding  SmS  and  2xSmS 
are  also  clearly  observed,  in  fact  2xSmS  is  detected  in  event  87183,  when  2xPmP  is 
not. 

Arrival  times  in  the  events  were  fitted  to  travel  time  curves  calculated  for  the 
P  and  S  phases  reflected  off  the  interfaces  of  the  two  cmstal  models  given  in  Table 
I.  They  were  obtained  for  the  E-W  trending  Flora-Asnes  profile  by  Sellevoll  and 
Warrick  (1971),  and  the  N-S  trending  Profile  3-4  by  Kanestrom  and  Haugland 
(1971).  The  locations  of  the  profiles  are  shown  in  Figure  7.  The  fitting  was  done 
by  adjusting  source  distances  and  origin  times  until  an  acceptable  match  in  travel 
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time  was  obtained.  Travel  time  curves  for  the  Flora-Asnes  profile  and  for  Profile 
3-4  are  plotted  with  the  record  sections  in  Figures  8  and  9,  respectively.  Distance 
adjustments  were  within  10  km  for  all  events  except  85313,  which  moved  from  a 
distance  of  251  km  to  287  km,  for  the  Flora-Asnes  profile,  and  305  km  for  Profile 
3-4.  An  acceptable  fit  was  obtained  for  both  models,  but  the  slopes  of  the  travel 
time  curves  for  Profile  3-4  better  match  the  observed  phase  velocities  of  the  seismo¬ 
grams;  the  slopes  on  the  PmP  and  SmS  curves  approach  the  limits  7.2  km/s  and  4.1 
km/s,  respectively.  Slopes  on  2xPmP  and  2xSmS,  however,  are  somewhat  lower 
than  the  observed  phase  velocities.  Slopes  on  the  Flora-Asnes  travel  time  curves 
are  generally  lower  than  the  observed  phase  velocities.  Travel  time  curves  for  P- 
to-S-converted  phases  are  also  included  in  Figure  8.  However,  they  do  not  correlate 
with  any  detectable  arrivals  in  the  composite  seismograms  indicating  that  the 
discontinuities  are  more  gradual  than  in  the  model. 

N-S  Section 

This  record  section  has  four  events,  which  after  relocation,  range  in  distance 
and  azimuth  from  73  to  170  km  and  157  to  202  degrees,  respectively;  the  two 
farthest  events  have  been  high-pass  filtered  before  stacking  (see  Fig.  10).  All  events 
have  a  clear  Pg  arrival,  with  a  phase  velocity  of  approximately  6.2  km/s,  but  a  less 
prominent  Sg.  PmP  is  very  faint  in  the  two  closest  events  (around  80  km),  which  is 
near  the  critical  distance,  but  becomes  clearer  with  increasing  distance.  Pg  is  the 
first  arrival  in  all  events,  but  PmP  has  started  to  interfere  with  it  in  event  88075,  at 
170  km  distance,  so  the  crossover  from  a  Pg  to  a  PmP  first  arrival  must  occur  just 
after  that  distance.  SmS  is  the  largest  arrival  in  all  but  one  event;  Event  88237  has 
an  Rg  wave  of  larger  amplitude  arriving  just  before  SmS.  A  relatively  smaller  Rg 
wave  arrives  before  SmS  in  event  88099.  Neither  Pn  nor  Sn  are  found  in  event 
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88154.  Event  88075,  on  the  other  hand,  has  a  faint  high-frequency  arrival,  which 
can  be  identified  with  Pn.  With  a  phase  velocity  around  9  km/s,  it  arrives  0.6 
seconds  before  PmP. 

Arrival  times  for  this  record  section  were  also  fitted  to  travel  time  curves  for 
the  two  crustal  models.  Both  models  can  match  the  observed  arrival  times  in  the 
records,  except  for  the  Moho  reflections  in  event  88099,  which  require  a  34-km- 
thick  crust.  The  slopes  of  the  travel  time  curves  for  the  slower  Flora-Asnes  model, 
however,  better  match  the  observed  phase  velocities.  The  travel  time  curves  for 
that  model  are  overlaid  on  the  record  section  in  Figure  10  The  crossover  from  a  Pg 
to  a  PmP  first  arrival  occurs  at  200  km  distance,  which  is  in  agreement  with  the 
phase  velocities  obtained  for  event  88075,  whereas  the  crossover  occurs  too  early  in 
the  Profile  3-4  model.  P-to-S-conversions  are  not  detected. 

The  effects  of  heterogeneous  structure  on  the  wave  propagation  along  the 
eastern  rim  of  the  Graben  are  frequency  dependent:  Waves  with  frequencies  less 
than  about  10  Hz  do  not  appear  to  be  affected  by  irregular  structure,  while  higher 
frequencies  tend  to  have  a  broad  power  maximum  in  the  wavenumber  plane,  with 
the  peak  stretched  toward  the  Graben.  Because  of  differences  in  frequency  content, 
the  four  events  are  affected  to  a  varying  degree.  The  average  frequency  spectra  of 
the  events  are  shown  in  Figure  11.  Due  to  the  low  signal-to-noise  ratio  in  events 
88154  and  88075  (approx.  4  and  2,  respectively)  their  spectra  peak  below  1  Hz. 
Above  that  peak  most  of  their  energy  is  carried  with  4.-5.  Hz  and  6.-8.  Hz  waves, 
respectively.  Most  of  the  energy  of  event  88099  is  concentrated  between  5  and  10 
Hz;  event  88237,  apart  from  the  1.1  and  1.6  Hz  Rg  waves,  has  most  of  its  energy 
carried  at  frequencies  between  8  and  14  Hz.  Due  to  their  deficiencies  in  high  fre¬ 
quencies,  events  88099,  88154,  ^d  88075  are  not  affected  by  heterogeneous 
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structure.  However,  event  88154,  due  to  its  closeness  to  the  Graben  area  may  exhi¬ 
bit  reflections  from  within  the  Graben;  the  two  small  phases  detected  between  the  P 
and  the  S  phases  have  backazimuths  of  213,  whereas  all  other  phases  have  backa- 
zimuths  around  200.  Event  88237,  due  to  its  high  frequencies,  is  greatly  affected  by 
irregular  crustal  structure.  Its  location,  however,  is  farthest  from  the  Graben.  P 
waves  with  peak  frequencies  above  10  Hz  generally  have  more  than  one  power 
maximum  in  the  wavenumber  plane,  at  backazimuths  around  130,  155,  and  170, 
whereas  lower  frequency  P  and  S  waves  have  only  one  at  backazimuths  between 
155  and  160.  When  relocating  the  event,  a  backazimuth  of  157  was  chosen. 

Discussion 

The  relatively  simple  structure  of  the  Baltic  Shield  is  the  dominant  cause  for 
the  clear  separation  and  sharpness  of  the  observed  crustal  arrivals  in  the  events  of 
both  record  sections  and  accounts  for  the  successful  modeling  of  the  arrivals  with 
homogeneous  plane  layered  models.  A  comparison  of  the  two,  however,  reveals 
some  differences  in  crustal  structure:  The  NW-SE  trending  profile  that  reaches  into 
the  Caledonian  mountains  appears  to  have  higher  crustal  velocities  than  the  N-S 
profile  that  runs  along  the  eastern  edge  of  the  Oslo  Graben.  Reflections  off  the 
Moho  are  sharper  near  the  critical  distance  and  have  higher  phase  velocities  in  the 
Caledonides  than  do  those  east  of  the  Graben,  indicating  higher  average  velocities 
in  the  crust  and  a  sharper  velocity  contrast  at  the  crust  mantle  interface.  Average 
crustal  thickness  appears  to  be  similar  for  both  profiles,  approximately  37  km,  but 
near  the  Graben  the  crust  thins  to  about  34  km.  This  is  in  agreement  with  previous 
studies,  which  have  shown  thinning  of  the  crust,  by  several  kilometers  in  the  Gra¬ 
ben  area  (Kanestrom,  1973,  Kanestrom  and  Haugland,  1971,  Cassel  et  al.,  1983). 
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The  high  phase  velcxities  obtained  for  Pn  also  indicate  a  dip  in  the  Moho  near  the 
array  of  approximately  10  degrees,  depending  on  velocity  model.  High  Pn  veloci¬ 
ties  under  the  NORSAR  array  have  also  been  reported  previously  by  Kanestrom  and 
Haugland  (1971)  and  Kanestrom  (1973).  Indeed,  velocities  in  a  NNW  direction 
were  found  to  be  10  km/s. 

Another  noticeable  difference  between  the  two  profiles  is  the  scattering  of 
waves  with  frequencies  above  10  Hz  in  the  N-S  profile.  This  scattering  is  not 
observed  in  the  NE-SW  profile.  To  illustrate  this  point.  Figure  12a  shows  the  dou¬ 
ble  power  peak  in  the  wavenumber  plane,  observed  in  the  Pg  phase  of  event  88237 
for  the  peak  frequency  of  10.8  Hz.  A  secondary  frequency  peak  in  the  same  time 
window,  but  of  lower  power  is  at  4.8  Hz.  This  lower  frequency  has  no  sign  of 
scattering,  as  shown  in  Figure  12b.  For  a  comparison  with  the  Caledonian  profile. 
Figure  12c  shows  the  power  spectrum  for  the  peak  frequency  in  the  Pg  window  of 
event  88267.  Its  amplitude  in  this  event  is  about  twice  the  amplitude  in  event 
88237.  At  a  frequency  of  13.7  Hz,  there  is  much  less  scattering  and  no  apptu-ent 
multipathing  as  in  the  southern  profile. 

For  a  comparison  of  the  coherency  in  the  events.  Figure  13  shows  plots  of  the 
envelope  stacks  (dashed  lines),  which  estimate  the  total  amplitude  in  each  event, 
and  the  incoherent  amplitude  (solid  lines)  for  four  of  the  events.  The  difference 
between  the  two  represents  the  coherent  amplitude.  Between  the  major  phase 
arrivals,  the  incoherent  amplitude  is  about  70  percent  of  the  total  ampliuide,  but  at 
the  arrivals  of  the  large  crustal  phases,  it  drops  down  to  about  20  percent  for  the  P 
phases,  somewhat  higher  for  the  S  phases.  The  incoherency  is  also  higher  for  the 
most  distant  event,  88286,  indicating  more  scattering  due  to  the  longer  travel  path. 
Event  88237,  which  showed  evidence  of  multipathing  has  high  incoherency 
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throughout  all  the  crustal  reflections  due  to  scattering  in  the  high-frequency  range, 
and  to  the  stacking,  which  favored  the  157  degree  backazimuth.  The  only  phase 
with  high  coherency  in  that  event  is  the  low-frequency  Rg  wave  at  about  20  sec.  on 
the  plot.  Event  88159  also  has  high  coherency  during  the  arrival  of  Rg.  Therefore, 
as  expected,  incoherency  increases  with  source  distance  and  frequency,  and  combin¬ 
ing  slant-stacks  to  approximate  a  seismogram  is  a  permitted  procedure,  since 
incoherency  during  the  arrivals  of  the  major  phases  is  low. 

Although  plane  layered  constant  velocity  models  can  account  for  all  major 
arrivals  in  the  seismograms,  data  from  events  at  a  distances  <  50  km  require  velo¬ 
city  gradients  near  the  surface.  This  is  in  accordance  with  the  refraction  study  of 
Mykkeltveit  (1980)  of  a  profile  north  of  the  Caledonian  cross  section  presented 
here.  A  future  improvement  of  the  structure  model  will  include  this  effect. 

The  events  presented  here  clearly  show  that  the  Pg  and  Lg  wave  trains  on 
regional  seismograms  have  separate  arrivals  within  them  with  distinctive  phase 
velocities,  which  constrain  their  travel  path,  and  at  the  same  time  the  crustal  struc¬ 
ture.  These  phases  can  be  valuable  in  modeling  structure,  as  well  as  in  event  loca¬ 
tion  and  source  depth  determination.  Describing  these  phases  only  in  terms  of  Pg 
and  Lg  is  therefore  an  oversimplification  and  leaves  valuable  information  contained 
in  them  unused. 
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TABLE  1. 


Flora-Asncs  Crustal  Model  (Sellevoll  and  Warrick,  1971).  • 

Vp  {km  sec  Vs  {km  sec  ') 

Density  {g  cm  ^) 

Thickness  (km) 

5.20  3.00 

2.60 

1 

6.00  3.46 

2.80 

14 

6.51  3.76 

3.00 

22 

8.05  4.65 

3.30 

Half-space 

Profile  3-4  Crustal  Model  (Kanestrom  and  Haugland,  1971).  * •• 

Vp  {km  sec'*)  Vs  {km  sec  *) 

Density  {g  cm  ^) 

Thickness  (km) 

5.20  3.00 

2.60 

1  •* 

6.20  3.58 

-A 

16 

6.55  3.78 

00 

10 

7.20  4.16 

3.13 

10 

8.20  4.73 

3.35 

Half-space 

*  S  wave  vcicKilics  arc  calculated  assuming  a  Poissons  solid. 

••  Surface  layer  added  to  the  model. 
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Figure  1 :  Map  showing  event  locations  (triangles)  and  location  of  NORESS 
(square).  Shaded  area  represents  the  Oslo  Graben. 
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Figure  2:  Composite  seismogram  for  event  A-6  and  time  sections  from  the  array 
beams  used  to  make  up  the  composite.  Bottom  trace  shows  the  whole 
2xSmS  beam. 


Figure  3:  Composite  seismograms  for  events  A- la  and  A- lb  plotted  at  the  source 
distance  determined  by  NORESS.  The  traces  are  normalized  to  their 
maximum  amplitude.  The  stacking  windows  are  indicated  by  dashed 
lines,  and  the  stacking  velocities  are  shown  on  the  left  margin  of 
each  window.  Arrows  indicate  arrival  times  of  phases  as  listed  in 
the  NORESS  bulletin.  The  thick  arrows  represent  the  phases  picked 
by  NORESS  as  P  and  Lg  and  are  used  in  the  automatic  location. 
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Figure  4:  Composite  seismograms  for  events  A-2a,  A-2b,  A-2c  and  A-2d  plotted  at 
their  source  distance  from  NORESS.  Same  scheme  as  Figure  3,  except 
for  arrival  times  in  events  A-2a,  A-2b  and  A-2c,  which  were  recorded 
during  a  period  of  unknown  time  drift  of  <  2  sec,  occurring  in  the 
on-line  data.  The  arrow  bars  indicate  the  uncertainty  in  the  arrival 
times  and  the  thick  bars  represent  NORESS  automatic  time  picks. 


-19- 


-20- 


Source 


Receiver 


Moho 


Figure  6:  A  schematic  explaining  the  notation  and  travel  paths  of  the  major 
crustal  and  upper  mantle  phases.  P  waves  converted  to  S  waves  at  the 
crustal  and  mantle  interfaces  will  be  denoted  by  Pg-S  and  PmS, 
respectively. 
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R«locat«d  Ev*nta 


Figure  7:  Map  showing  relocated  events  (solid  triangles)  and  location  of  NORESS 
(solid  square).  Hatched  area  represents  the  Oslo  Graben.  The  NW-SE 
line  represents  the  Flora-Asnes  refraction  line,  and  the  N-S  line 
represents  the  Profile  3-4  refraction  line. 
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Section  A  —  Flora-Asnas  Modal 


Figure  8:  Composite  seismograms  for  the  events  from  the  Caledonian  section 
plotted  as  a  function  of  distance,  reduced  with  a  velocity  of  7  km/s 
and  normalized  to  maximum  amplitude.  Dashed  lines  above  seismograms 
indicate  where  the  slant-stack  sections  have  been  pieced  together. 

Stacking  velocity  is  shown  on  the  left  margin  of  each  interval. 

Travel  time  curves  for  the  Flora-Asnes  model  are  overlaid  on  the 
plot.  They  represent:  (— )  PmP  and  SmS;  Pg  and  Sg; 

^...)  Pg-S  and  PmS;  2xPmP  and  2xSmS;  (-- )  2xPg.  Slopes 

representing  velocities  on  the  curves  can  be  read  from  the  velocity 
template  to  the  right  of  center  on  the  plot. 
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Figure  9:  Same  record  section  as  in  Figure  8,  but  with  the  travel  time  curves 
for  the  Profile  3-4  model  overlaid.  Reflections  off  the  two  inter¬ 
faces  in  the  crust  are  referred  to  as  Pgl  and  Pg2,  respectively. 

A  similar  notation  is  used  for  S-waves.  Travel  time  curves 
represent:  ^  Pgl  and  Sgl;  Pg2  and  Sg2.  Same  notation  for 
Moho  reflections  as  in  Figure  8. 
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Diataaoc  (km) 


Figure  10:  Composite  seismogram  record  section  from  the  N-S  profile,  with  the 
travel  time  curves  for  the  Flora- Asnes  profile  overlaid.  Same 
scheme  as  in  Figure  8  except  travel  time  curves  for  the  double  Moho 
reflections,  2xPmP  and  2xSmS  are  not  shown. 
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Figure  1 1 :  Average  amplitude  spectra  of  the  records  from  the  four  events  in  the 
N-S  record  section,  each  calculated  for  a  30  sec  time  window  con¬ 
taining  tho  whole  waveform. 


NBFK  -  lingle  window 
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NBFK  -  single  window 
Estimated  bearing;  155. 

Estimated  velocity:  6.41 

Analysis  frequency;  4.800 
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Scaling  type:  LINEAR 
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NBFK  -  single  window 
Estimated  bearing:  303. 

Estimated  velocity:  7.24 

Analysis  frequency;  13.70 
Max.  horiz.  wavenumber;  2.5000 
Scaling  type:  UNEAR 
Max.  power  value:  O.BASGBe^-OB 


NORTH 


Figure  12:  Power  spectra  at  peak  frequencies  in  Pg  for  events  B-2  and  A-2d. 

a)  Event  B-2  at  maximum  peak  frequency. 

b)  Event  B-2  at  secondary  peak  frequency. 

c)  Event  A-2d  at  maximum  peak  frequency. 
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Figure  13:  Estimated  total  amplitude  (dashe 
for  events  A-5d,  A-3,  A-2c  and  B-2. 
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Abstract 


Broad- band,  digital,  teleseismic  P-wave  data  from  the  University  of  California,  Berkeley, 
stations  BKS,  MHC,  and  SAO  are  used  to  study  the  character  of  teleseismic  P  coda.  Velocity 
heterogeneity  in  structure  near  these  stations  gives  rise  to  a  variety  of  scattering  effects  which  include 
azimuthly  dependent  receiver  functions,  large  tangential  component  particle  motions,  and  high 
levels  of  P  coda.  P  codas  for  shallow  and  deep  telcseisms  are  analyzed  for  three  frequency  bands 
centered  at  1.5,  0.75,  and  0.375  hz.  Observed  coda  decay  rate  and  level  are  fit  with  an  energy  flux 
theory  to  parameterize  scattering  effects  near  the  source  and  at  the  receiver.  Coda  decay  rate  and 
level  for  deep  earthquakes  are  seen  to  be  stable  observables  between  stations,  implying  that  the 
scattered  coda  field  is  homogeneously  distributed  over  the  receiver  network  validating  a  major 
assumption  in  the  energy  flux  model.  Scattering  Q,  Q,,  is  not  particularly  well  constrained  by  the 
data  and  is  inferred  to  be  affected  by  focussing  and  defocussing  effects  of  the  direct  P-waves.  There 
is  a  pronounced  difference  in  coda  decay  rate  between  shallow  and  deep  earthquakes  for  1.5  hz 
coda  where  coda  decay  rate  for  deep  events  is  about  a  factor  of  two  greater  than  for  shallow  events. 
There  is  no  difference  in  coda  decay  rate  between  the  two  populations  of  events  in  the  lower  fre¬ 
quency  bands.  The  observation  at  1 .5  hz  is  consistent  with  a  simple  model  of  P-wave  coda  formation 
where  most  coda  forms  in  a  scattering  layer  near  the  surface  of  the  earth  in  regions  near  the  .source 
and  near  the  receiver.  The  homogeneity  of  coda  decay  at  0.75  and  0.375  hz  for  deep  and  shallow 
events  suggests  that  a  large,  near-receiver  scattering  mechanism,  such  as  surface  wave  formation, 
dominates  the  coda. 
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Introduction 

Scattering  of  teleseismic  body  waves  in  heterogeneous  structure  is  a  source  of  signal-induced 
noise  for  many  broad-band  wave  form  studies.  Efforts  to  deterministically  model  the  wave  form 
rely  on  specific  knowledge  of  earth  structure  that  can  be  represented  in  reasonably  simple  plane 
layered  or  vertically  inhomogeneous  models.  Most  of  the  effort  expended  in  detailed  source 
mechanism  studies  lies  in  determining  the  tradeoffs  between  inferred  source  parameters  and 
unknown,  but  bounded  estimates  of  earth  structure.  Signal  produced  by  waves  scattered  by  lateral 
heterogeneity  will  naturally  fold  back  into  erroneous  models  of  the  source  since  the  wave  propa¬ 
gation  will  be  assumed  to  be  known  through  synthetic  seismogram  calculations  for  simple  structure. 
The  scattered  signal  may  also  be  inadvertently  modeled  for  simple  plane  layered  earth  structure  in 
receiver  studies,  yielding  misleading  crustal  or  mantle  structures. 

A  principal  characteristic  of  scattering  is  the  production  of  an  extended  coda  after  major 
arrivals.  Many  studies  have  been  performed  to  study  coda  particularly  for  cases  of  local  S  wave 
propagation  (see  Herraiz  and  Espinosa,  1987,  for  a  review).  The  obvious  effect  of  scattering  in 
teleseismic  P-waves  is  to  produce  an  extended  coda  with  a  characteristic  exponential  decay.  The 
extended  coda  becomes  an  important  contaminant  if  there  are  later  arriving  phases  of  interest  such 
as  small  Ps  conversions  from  crustal  interfaces.  However,  the  exponential  decay  implies  that 
scattered  waves  will  be  largest  at  the  P-wave  arrival  itself.  Becau.se  these  scattered  waves  arrive 
netu  Fermat’s  minimum  time,  they  can  also  give  rise  to  focussing  or  defocussing  effects,  depending 
on  the  specific  structural  geometry,  which  can  severely  distort  the  incident  P  wave  front. 

A  previous  study  (Langston,  1989)  examined  scattering  seen  in  teleseismic  P  wave  receiver 
functions  from  deep  events  recorded  at  PAS  (Pa.sadena,  California)  and  SCP  (State  College, 
Pennsylvania).  Scattering  from  heterogeneous  structure  at  PAS  was  evident  through  observations 
of  large,  off-azimuth  body  wave  conversions  and  high  levels  of  coda  over  the  three  components  of 
motion.  Comparison  of  codas  observed  at  PAS  and  SCP  showed  that  high  levels  of  coda  correlate 
with  slow  coda  decay  rate.  These  observations  were  modeled  using  an  energy  flux  model  for 
scattered  body  waves  in  a  heterogeneous  layer  over  an  isotropic  halfspace.  Simulations  of  body 


wave  propagation  in  random  one-dimensional  (ID)  structure  showed  that  observed  codas  were  too 
large  to  be  explained  by  plausible  ID  structures.  However,  the  ID  simulations  did  show  that  as 
the  severity  of  scattering  increased  with  heterogeneity,  coda  decay  rate  decreased  as  coda  level, 
relative  to  the  direct  P  wave,  increased,  in  accord  with  the  data.  If  interpreted  with  a  single  .scattering 
model,  the  slower  coda  decay  rate  would  imply  a  higher  scattering  Q  or  less  scattering  attenuation. 
The  apparent  paradox  was  resolved  with  the  energy  flux  model  where  scattering  Q,  as  usually 
defined,  was  decoupled  from  the  process  which  controls  coda  decay.  Coda  decay  rate  is  controlled 
by  radiative  diffusion  of  the  coda  field  from  the  scattering  layer  into  the  isou-opic  halfspace. 

Although  the  model  allows  a  simple  parameterization  of  teleseismic  P  coda,  it  primarily 
introduces  some  interesting  issues  in  both  the  interpretation  of  coda  and  its  mechanism  of  formation. 
Observed  coda  levels  in  the  receiver  function  data  are  several  times  higher  than  expected  from  body 
wave  scattering  and  suggest  that  body-to-surface  wave  conversion  is  important.  Blakeslee  and 
.Malin  (1990)  studied  small  earthquakes  at  Parkfield,  California,  with  a  borehole  array  and  saw  that 
coda  complexity  increased  near  the  surface.  They  suggested  that  body  wave  scattering  dominated 
coda  away  from  the  free  surface  and  that  body  wave-to-surface  wave  conversions  dominated  coda 
at  the  surface.  This  mechanism  of  coda  production  will  be  frequency  dependent.  Considering  the 
reciprocal  problem  of  a  source  in  the  scattering  layer,  it  is  also  likely  that  near-  source  coda  production 
follows  similar  mechanisms  where  surface  wave-to-body  wave  interactions  occur  and  may  be  source 
depth  dependent  (McLaughlin  et  al,  1985).  An  expectation  of  this  simple  view  of  coda  formation 
is  that  shallow  earthquakes,  as  seen  teleseismically,  should  display  higher  levels  of  coda  with  slower 
decay  rate  than  seen  for  deep  earthquakes.  Dainty  et  al  ( 1974)  suggested  this  kind  of  model  as  an 
explanation  of  lunar  seismogram  coda  using  an  energy  diffusion  approach  for  scattering  in  the  lunar 
regolith.  Thus,  it  might  be  expected  that  the  apparent  coda  Q  for  shallow  events  should  generally 
be  greater  than  that  measured  from  deep  events  if  scattering  is  confined  by  heterogeneity  in  the 
lithosphere. 


A  major  assumption  in  the  energy  flux  coda  model,  as  presented  by  1  rankel  and  Wennerberg 
(1987)  for  a  wholespace  and  modified  for  a  scattering  layer  over  a  homogeneous  halfspace  by 
Langston  (1989a,  1989b)  and  Kom(1990),  is  that  the  scattered  coda  field  arranges  itself  homoge¬ 
neously  behind  the  incident  wave  front.  Frankel  and  Wennerberg  (1987)  developed  the  model  by 
appealing  to  finite  difference  simulations  of  wave  propagation  in  scattering  media.  Although  they 
show  that  this  assumption  is  very  good  in  the  simulations,  empirical  verification  with  the  data  is 
desirable.  Aki  and  Chouet  (1975)  suggest  in  their  study  of  local  S  wave  codas  that  this  is  true  since 
the  apparent  mean  free  path  of  scattered  coda  waves  is  large  compared  to  the  travel  path  of  the 
direct  S  wave.  A  purpose  of  this  study  is  to  use  data  from  a  regional  network  to  see  if  P  coda  displays 
similar  characteristics  of  decay  and  level  over  the  region,  substantiating  this  assumption  of  coda 
homogeneity. 

In  this  paper,  we  use  data  collected  from  the  University  of  California  Berkeley  seismic  network 
U)  examine  coda  from  a  suite  of  tele.sei.smic  events  at  several  receivers.  The  data  were  collected  as 
part  of  a  receiver  function  study  to  determine  gross  crustal  and  upper  mantle  structure  in  Central 
California.  This  aspect  of  the  study  will  be  reported  elsewhere.  We  wish  to  primarily  investigate 
the  stability  of  coda  decay  and  level  estimates  by  examining  coda  from  the  same  events  recorded 
at  different  receivers.  The  results  will  show  that,  like  many  coda  determinations  from  local  S  waves. 
c(xla  decay  rate  is  a  generally  stable  measure  for  individual  events.  However,  coda  decay  rate  also 
depends  on  source  depth;  being  slower,  in  general,  for  shallow  sources.  This  supports  a  simple 
model  of  scattering  in  the  earth  suggested  by  Dainty  and  Harris  (1989)  where  scattering  occurs 
primarily  at  the  source  and  at  the  receiver  with  little  scattering  occurring  along  the  mantle  path  of 
the  seismic  ray.  The  results  will  also  show  that  inferring  scattering  Q  from  the  relative  coda  level 
is  generally  unstable  and  is  probably  contaminated  by  direct  wave  focussing  effects. 


Berkeley  Network  Data 


Ammon  (1985)  analyzed  deep  eanhquake  P-wave  forms  recorded  on  the  analog  system  at 
BKS  to  recover  the  receiver  function  for  structure  under  the  station.  He  found  that  the  radial 
component  receiver  functions  varied  substantially  with  azimuth  and  contained  many  large  arrivals 
compared  to  direct  P.  Particle  motions  were  also  complex  for  later  arrivals  giving  rise  to  tangential 
motions  relative  to  the  P-wave  great  circle  path.  An  original  purpose  of  the  present  study  was  to 
obtain  high  quality  digital  data  from  Berkeley’s  new  broad-band  digital  system  to  examine  these 
scattering  effects  in  more  detail  and  to  extend  structure  studies  to  other  receivers  in  the  region. 
Figure  1  shows  the  location  of  BKS  (Berkeley),  MHC  (Mount  Hamilton)  and  SAO  (San  Andreas 
Observatory)  stations  which  have  broad-band,  three  component,  digital  instrumentation.  Bolt  et  al 
( 1988)  and  Collins  and  LIhrhammer  (1988)  describe  the  digital  recording  system  and  responses  for 
these  stations.  Responses  are  flat  to  ground  velocity  between  0.05  and  5  hz.  Earthquake  events 
were  chosen  for  study  by  examining  the  catalog  of  available  waveforms  at  the  Berkeley  Seismic 
Station  and  choosing  those  with  go<xl  signal  to  noise  characteristics.  Usually  this  meant  finding 
teleseisms  with  body  wave  magnitude  greater  than  5.5.  Table  1  contains  event  parameters  and 
Table  2  the  stations  recording  each  event. 

Figure  2  shows  the  three  components  of  motion  near  the  P-wave  arrival  for  the  10/08/88 
Tonga  Islands  event  recorded  at  BKS.  MHCand  BKS  have  the  same  general  receiver  characteristics 
showing  azimuthaliy  varying  radial  receiver  functions  and  a  tangential  component  about  20%  of 
the  size  of  the  direct  P  wave  arrival.  The  arrow  (Figure  2)  points  to  an  anomalously  large,  later- 
arriving  receiver  phase  on  the  radial  component  which  is  substantially  larger  than  direct  radial  P. 
This  arrival  is  not  seen  in  other  azimuths  implying  that  plane  layered  earth  structure  may  not  be  a 
good  approximation  for  structure  in  the  area  (e.g.,  Langston,  1979). 

Figure  3  shows  P-wave  form  data  recorded  at  SAO  from  the  8/14/88  Sea  of  Okhotsk  event. 
SAO  structure  is  inferred  to  be  radically  heterogeneous  near  the  surface  since  the  tangential  com¬ 
ponent  of  motion  for  first  arriving  P  is  as  large  as  the  vertical  component.  This  effect  changes  with 


azimuth  and  frequency  and  will  be  described  in  the  companion  study.  Since  SAO  is  located  at  the 
trace  of  the  San  Andreas  fault,  the  extremely  large  near-surface  effect  on  panicle  motions  must  be 
due  to  velocity  structure  within  the  fault  zone  (e.g.,  Thurber,  1983).  Of  interest  to  this  study  on 
ctxla  is  the  simple  observation  that  the  direct  P  wave  is  severely  affected  by  lateral  heterogeneity 
and  may  be  focussed  or  defocussed.  The  amplitude  of  direct  P,  relative  to  later  arriving  coda  waves, 
is  important  when  inverting  for  coda  parameters. 


Theory  and  Data  Analysis 

Details  of  the  development  of  the  energy  flux  coda  model  are  found  in  Langston  (1989a). 
Conceptially,  the  model  is  based  on  the  assumption  of  a  vertically  incident  P  plane  wave  which 
propagates  upward  from  a  homogeneous  halfspace  into  a  layer  containing  scatterers.  The  P-wave 
is  assumed  to  be  impulsive  so  that  its  duration  is  much  less  than  the  traveltime  through  the  layer. 
Consideration  of  conservation  of  energy  scattered  from  the  P  wave  into  coda  gives  the  level  ol  coda 
amplitude.  Ignoring  the  displacement  of  the  direct  P  wave,  coda  amplitude  at  the  free  surt'ace  of 
the  layer  is  given  by: 


where,  w  is  circular  frequency,  is  the  average  travel  time  of  the  direct  wave  through  the  layer, 
is  the  scattering  Q  of  the  layer  and  Q,  is  the  "radiative  "  Q  which  governs  diffusion  of  coda  waves 
out  of  the  layer  into  the  underlying  halfspace.  Given  ground  velocity,  v(t),  of  the  incident  P  wave 
pul.sc  centered  around  some  circular  frequency  w.  the  factor  /„  is  given  by 


(2) 


The  integration  bound  t,  is  the  P  wave  arrival  time  but  tj  must  be  estimated  from  the  data  and  must 
be  such  that  it  includes  all  of  the  incident  P  pulse.  Estimating  tj  is  difficult  for  emergent  and  long 
duration  signals. 

ITie  theory  for  equation  (1)  was  developed  assuming  a  scalar  wave  field,  such  as  pressure, 
with  scattering  confined  to  the  same  wave  type.  Equation  (1)  is  applied  to  three-component  data 
with  the  understanding  that  the  resulting  parameterization  results  in  empirical  measures  of  coda 
decay  and  level.  No  assertions  are  made  on  the  scattering  mechanisms  that  go  into  coda  fomiation. 
The  three-component  P-wave  data  for  an  individual  earthquake  are  prtKessed  in  the  following 
manner  since  coda  occurs  on  each  component.  First,  each  comjxinent  is  bandpass  filtered  using  a 
zero  phase  Butterworth  filter  centered  at  a  given  frequency.  The  factor /n  is  computed  by  choosing 
t|  and  tj  from  the  vertical  component  of  motion  since  the  P  pulse  is  most  clear  on  this  component. 
The  integral  is  computed  for  the  sum  of  the  squares  of  the  three  components  over  the  inferred  time 
interval.  This  is  the  value  of  Ip.  Next,  the  amplitude  envelope  of  each  component  is  found  by 
computing  the  modulus  of  the  analytic  signal.  The  three  envelopes  are  combined  by  summing  the 
square  of  each  amplitude  sample  point  by  point.  The  final  coda  envelope  is  the  square  root  of  the 
sum  and  is  normalized  by  dividing  by 

Once  a  coda  envelope  is  obtained  for  an  event,  equation  (1)  is  fit  to  the  envelope  using  a 
least-squares  algorithm  on  the  logarithm  of  the  coda  envelope.  The  logarithm  is  used  since  the 
inversion  will  be  linear  for  Qr\  The  zero  intercept  of  the  log  coda  envelope  is  related  to  Q,"'  in  a 
complicated  way,  but  the  least-squares  inversion  is  performed  iteratively  until  convergence. 
Synthetic  tests  showed  that  there  was  little  starting  model  dependence  in  obtaining  the  0“' 
parameters  and  convergence  occurred  in  all  cases  within  five  iterations. 

Filtering  was  performed  at  three  different  center  frequencies  based  on  trial  processing  and 
examining  noise  levels  in  the  teleseismic  data.  These  center  frequencies  were  at  1.5  (comer  fre¬ 
quencies  of  1  and  2  hz),  0.75  (comer  frequencies  of  0.3  and  1  hz)  and  0.375  hz  (comer  frequencies 
at  0.25  and  0.5  hz).  A  three  pole,  phaseless,  Butterworth  filter  was  used  in  each  case. 


Zero  lapse  time  for  the  coda  was  taken  at  the  peak  of  the  P  wave  arrival.  This  must  be  known 
to  apply  equation  (1).  The  Q  parameters  were  determined  from  the  coda  envelope  for  each  event 
by  inversion  of  80  seconds  of  data  starting  20  seconds  after  the  P  wave  peak  to  avoid  the  direct  P 
wave.  Tables  3,  4  and  5  contain  the  values  for  each  event  and  station,  t^  was  taken  to  be  10  sec 
using  an  average  P  wave  velocity  of  6  km/s  for  a  scattering  layer  60  km  thick. 

Discussion 

The  average  standard  deviation  on  inferred  Q  parameters  was  about  20%.  However,  it  was 
found  that  coda  decay  rate  measured  by  Q,  depended  on  source  depth  and  coda  level  inferred  from 
Os  depended  on  recording  station.  Figure  4  shows  correlograms  of  Q,  and  Qj  for  the  three  stations. 
The  plots  are  for  coda  parameters  taken  from  events  which  were  recorded  at  the  starion  pair  BKS 
and  MHC  and  at  the  pair  SAO  and  MHC.  Although  values  for  Q,  scatter  in  magnitude  among  the 
event  data,  the  correlograms  show  that  Q,  correlates  well  between  stations  for  the  same  event. 
Equation  (1)  shows  that  Q,  is  proportional  to  the  inverse  slope  of  the  log  coda  envelope.  Thus,  the 
giKxl  correlation  for  Q,  shows  that  ctxla  decay  rate  is  a  stable  quantity  over  the  region  for  a  single 
event. 

The  scattering  Q,  Q„  is  less  well-behaved.  Figure  4  shows  that  the  correlation  between  stations 
IS  not  very  good  and  that  Q,  dciemiinaiions  can  easily  vary  by  factors  of  two  for  the  same  event. 
Qj  is  a  measure  of  coda  level  relative  to  the  direct  wave  and  therefore  depends  directly  on  how  well 
yJTp  is  determined.  A  major  assumption  in  the  coda  theory  is  that  the  direct  wave  only  loses  energy 
in  scattering  and  is  not  affected  by  any  other  wave  propagation  effect.  Thus,  the  reliance  of  coda 
level  on  the  direct  wave  implicitly  assumes  that  the  direct  wave  is  not  focussed  or  defocussed  by 
the  scattering  medium.  This  is  probably  the  weakest  assumption  contained  in  any  coda  theory  since 
the  amplitude  of  the  input  wave  to  the  scattering  medium  is  never  a  known  quantity.  Figure  5 
demonstrates  that  differential  focussing  of  the  direct  wave  occurs  over  the  network.  Shown  are  the 
P-wave  vertical  components  for  the  9/07/88  Honshu  event.  Note  the  factor  of  three  difference  in 


amplitude  between  BKS  and  SAO  for  the  direct  P-wave  yet  the  coda  levels  remain  comparable. 
The  empirical  evidence  of  Figures  4  and  5  suggests  that  determination  of  coda  level  and  Q,  is  a 
noisy  process  in  the  analysis  of  teleseismic  P  codas  at  single  receivers  and  that  it  may  be  better  to 
forgo  normalization  by  the  direct  wave  amplitude.  Constancy  seen  in  ctxla  levels  across  the  network 
also  strongly  supports  the  assumption  of  general  homogeneity  of  the  scattered  coda  field,  justifying 
this  assumption  in  the  energy  flux  model. 

Figure  6  shows  a  compilation  of  all  inferred  Q,  values  plotted  versus  event  source  depth.  High 
Q,  values  imply  that  slope  is  small  and  that  coda  decays  slowly.  In  the  single  scattering  model  (Aki 
and  Chouet,  1975),  Q,  is  interpreted  in  terms  of  scattering  attenuation;  high  slopes  mean  greater 
attenuation.  However,  Langstonf  1989a)  suggested  that  slow  coda  decay  rate  really  implies  that 
scattering  effects  are  more  severe  than  .seen  in  teleseismic  P  codas  which  decay  more  quickly. 
Empirically,  slow  P-coda  decay  rates  are  associated  with  high  coda  levels  relative  to  the  direct  P 
wave  arrival. 

The  data  (Figure  6)  show  that  coda  decay  rate  can  be  source  depth  dependent.  Results  for 
the  1.5  hz  bandpass  show  that  coda  decay  rate  is  about  half  as  fast  for  the  events  above  100  km 
depth  compared  to  deeper  events.  The  data  for  the  other  two  bandpasses  do  not  show  the  same 
effect,  however.  It  appears  that  coda  decay  rate  is  comparable  for  all  events,  regardless  of  source 
depth,  for  the  0.75  hz  and  0.375  hz  bandpasses.  Oxla  data  for  all  deep  events  at  all  stations  were 
stacked  to  get  the  best  estimate  of  Q,  and  Q,  for  the  central  California  region.  This,  of  course, 
assumes  that  coda  associated  with  scattering  near  the  source  for  these  deep  events  is  a  minor  con¬ 
tributor  to  the  observed  coda,  an  assumption  that  is  discussed  below.  The  coda  and  model  fits  are 
shown  in  Figure  7.  Note  the  linearity  of  the  coda  slopes  and  how  coda  slope  does  not  change  much 
between  the  three  frequency  bands.  Coda  levels  and  slopes  are  comparable  to  those  found  for  PAS 
station  in  Southern  California  for  teleseismic  deep  events  (Langston.  1989). 

A  simple  model  for  coda  formation  is  that  most  scattering  occurs  in  the  lithosphere  with  little 
or  no  scattering  occurring  along  the  body  wave  ray  path  throughout  most  of  the  deeper  mantle. 
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Thus,  coda  is  produced  at  the  receiver  from  the  random  scattering  of  waves  from  lateral  hetero¬ 
geneities  and  probably  consist  of  body  waves  and  body-to-surface  wave  conversions  (e  g.  Blakeslee 
and  Maiin,  1990;  Langston,  1989a).  Coda  waves  are  also  produced  at  the  sodct,  for  sources  within 
the  lithosphere,  and  because  of  the  teleseismic  distances  involved,  travel  to  the  receiver  as  P  w  aves 
along  paths  very  close  to  the  direct  P  wave  arrival.  Although  near-so-T-^  scattering  can  involve 
complex  effects,  including  for  example  Rayleigh  wave-to-P  wave  conversions,  the  constraint  on 
coda  arrival  times  being  comparable  to  the  direct  wave  guarantees  that  source  coda  simply  be  a  pan 
of  the  effective  source  time  function  when  it  arrives  at  a  teleseismic  receiver. 

The  relative  dominance  of  source  and  receiver  coda  can  be  assessed  by  the  following  simple 
analysis.  Coda  at  the  receiver  can  be  thought  of  as  the  result  of  convolving  the  response  of  the 
receiver  function,  r(t),  with  the  effective  source  function,  s(t).  The  deep  event  data  suggest  that 
coda  slope  and  approximate  levels  do  not  change  much  with  frequency.  Consider  the  following 
functional  form  for  the  source  and  receiver  functions: 

.v(f)  =  5(t)-i-C,(/)c  (3) 

r(t)^d(t)  +  C,{t)e  (4) 

where  C,(t)  and  C!,(t)  are  random  time  series  given  by 

C,(r)  =  C,  S  a  dd  -JAt)  (5) 

/  =0 

C^^)  =  C,^  b^bd-kAi)  (6) 

*  =n 

and  C,  are  source  and  receiver  coda  scaling  factors  controlling  coda  level  and  y,  and  y,  are  the 

■source  and  receiver  ctxla  decay  rate  parameters  for  the  source  and  receiver  ccxlas.  respectively.  The 
coefficients  a^  and  bj  of  equations  (5)  and  (6)  arc  Gaussian  distributed  variates  between  -  I  and  1. 
with  zero  mean,  sampled  every  At.  The  Dirac  delta  functions,  6(t ),  in  equations  ( 3)  and  (4)  represent 


the  "deterministic"  response  of  the  direct  P  wave  from  the  source  and  at  the  receiver.  The  relative 
size  of  coda  is  controlled  by  the  scaling  factors  and  decay  rates.  Convolving  (3)  and  (4)  gives  the 
vertical  component  of  displacement,  w(t)  (see  Appendix) 
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w(r)  =  m  +  C,(t)e^‘'  +  CXt)e^'‘  +  C,C,«(r)— (7) 

Jr  Is 

for  different  decay  parameters  and 

w(r)  =  m  +  {C,(r)  +  C,(r)}e-^  +  C,C,u(r)re-^  (8) 

if  tl  e  two  decay  parameters  are  the  same  and  equal  to  y.  The  function  u(t)  is  also  a  Gaussian- 
distributed  random  time  series. 

According  to  (7),  if  one  of  the  decay  constants  is  substantially  less  than  the  other  (implying 
more  scattering  or  higher  levels  of  coda)  for  approximately  equal  initial  levels,  then  the  slope  of 
the  log  coda  envelope  will  be  largely  determined  by  this  slower  decay  rate.  If  the  decay  constants 
are  comparable  or  the  same  (equation  (8)),  then  the  coda  will  have  the  same  amplitude  decay  at 
long  lapse  times  but  at  shorter  times  may  be  dominated  by  the  additional  factor  of  i  in  the  last  term. 
In  the  case  of  the  deep  event  data  studied  here  (Figure  7)  y,  =  0.023  sec  '  on  average  for  the  three 
frequency  bands  (assuming  that  C,  =  0).  If  an  effective  source  function  from  a  shallow  event  has 
a  coda  decay  rate  comparable  to  this  receiver  coda  decay  rate,  then  the  factor  with  t  in  equation  (8) 
may  be  of  the  same  order  as  the  other  exponential  in  the  log  coda.  This  has  the  effect  of  flattening 
the  log  coda  envelope  function  biasing  the  slope  to  smaller  decay  constants. 

An  examination  of  the  shallow  event  coda  data  for  the  1.5  hz  bandpass  showed  both  simple 
linear  coda  decay  starting  at  the  P  arrival  for  some  events  and  more  complex  envelopes  mimicking 
the  behavior  implied  by  equation  (8)  for  other  events.  It  should  be  stated  that  near-source  rever¬ 
berations  and/or  complex  and  lengthy  source  time  histories  for  the  shallow  events  can  contribute 


to  increasing  apparent  coda  levels  at  relatively  short  lapse  times.  It  may  be  possible  that  these 
additional  factors  contribute  to  the  difference  in  coda  decay  rates  seen  between  the  deep  and  shallow 
events. 

However,  there  is  no  a  priori  reason  why  there  should  be  a  difference  in  Q,  parameters  between 
shallow  and  deep  events  at  1 .5  hz  and  not  for  the  other  frequency  bands  (Figure  6).  The  analysis 
technique  of  stacking  coda  from  all  three  components  of  motion  at  the  receiver  will  bias  results 
toward  the  receiver  since  each  horizontal  component  contributes  just  as  much  energy  to  the  receiver 
coda  as  the  vertical  component  and  is  fundamentally  a  result  of  scattering  near  the  receiver.  Thus, 
source-area  scattering  must  be  large  to  overcome  or  substantially  add  to  the  receiver  coda  signal. 
The  data  of  Figure  6  suggests  that  scattering  is  very  efficient  near  shallow  sources  at  1.5  hz.  Perhaps 
the  difference  between  frequency  bands  is  related  to  a  change  in  scattering  mechanism  such  as 
increased  Rayleigh-to-P  wave  conversion. 

That  some  of  the  data  show  a  source  depth  dependent  coda  decay  rate  suggests  that  source 
area  .scattering  is  important  in  forming  tele.seismic  P  coda.  However,  the  ambiguity  remains  con¬ 
cerning  the  relative  contribution  of  receiver  and  source  area  scattering.  Several  lines  of  evidence 
suggest  that  receiver  structure  heterogeneity  u.sually  dominates  the  formation  of  coda,  panicularly 
when  deep  event  data  are  considered.  First,  particle  motions  are  observed  to  be  complex  within 
the  coda  with  very  little  evidence  of  rectilinear  motions  characteristic  of  body  waves.  Tangential 
motions  almost  always  attain  similar  amplitude  levels  to  radial  motions  in  ctxla  after  teleseismic  P. 
This  remains  true  even  after  a  source  equalization  operation  where  the  vertical  component  is 
deconvolved  from  the  horizontals  in  receiver  function  studies.  When  data  from  the  same  event  are 
considered  at  several  stations,  as  in  Figure  5,  focussing  effects  are  confined  to  the  early,  impulsive 
parts  of  the  record  with  the  ctxla  level  being  unaffected  between  stations.  This  has  also  been  a 
common  observation  at  teleseismic  arrays  which  show  large  (order  of  magnitude)  direct  wave 
amplitude  fluctuations  for  short-period  P  waves  between  elements  of  the  array  (e.g.,  Bungum  et  al, 
1971 ).  If  significant  coda  were  contained  within  the  incident  P  wave  from  the  source  area,  it  would 
be  expected  that  any  ffx'iissing  effects  would  identically  affect  both  the  impulsive  part  of  the  P  wave 


as  well  as  the  source  area  coda  portion  of  the  incident  P.  Since  this  is  not  observed  for  the  deep 
events  and  because  of  obvious  tangential  components  indicative  of  near-receiver  scattering,  we 
ascribe  coda  decay  and  level  determined  from  the  deep  event  data  (Figure  7)  to  be  due  to  scattering 
near  the  three  receivers.  The  good  correlation  of  coda  decay  rate  and  absolute  coda  level  for  the 
same  deep  event  at  the  three  receivers  strongly  supports  the  assumption  of  a  homogeneously  dis¬ 
tributed  coda  field  required  by  the  energy  flux  model. 

Conclusions 

Coda  induced  by  scattering  of  teleseismic  body  waves  along  the  Hayward-Calaveras  fault 
system,  a  structurely  complex  area  of  central  California,  displays  characteristics  simihu-  to  coda 
seen  in  receiver  function  data  at  Pasadena,  California.  Coda  levels  are  relatively  high  compared  to 
the  direct  P  wave  and  the  decay  rate  is  slow.  The  use  of  broad  band  data  from  several  receivers  in 
the  Berkeley  network  allow  comparison  of  the  stability  of  coda  decay  rates  and  levels  over  a  fre¬ 
quency  band  of  2  to  0.25  hz.  Using  an  energy  flux  model  to  parameterize  coda  decay  rate  by  a 
radiative  Q  piirameter,  Q„  and  the  coda  level  by  the  scattering  Q,  Q,,  the  data  in  Central  California 
>how  that  coda  decay  rate  is  largely  frequency  independent  and  consistent  between  stations.  Ctxla 
levels  are  referenced  to  the  direct  P  wave  and  Q,  values  show  larger  amounts  of  scatter  between 
stations  due  to  foccussing  and  defocussing  of  the  direct  P  wave  by  heterogeneity  under  the  receivers. 
Cf'da  decay  rate  for  1.5  hz  bandpass  data  from  shallow  earthquake  source.,  is  about  twice  as  slow 
as  coda  decay  rate  from  deep  events.  This  effect  is  not  seen  at  lower  frequencies.  The  slower  coda 
decay  rate  implies  that  the  contribution  of  coda  due  to  source  area  scattering  increases  above  1  hz 
and  is  comparable  to  scattering  effects  at  the  receivers. 


Table  1 

Event  Parameters 


Event 

Date 

Origin  Time 

Lat.  (“) 

Long.(“) 

Depth 

(km) 

mb 

Region 

1 

02/10/87 

00:59:28.5 

19.49S 

117.45E 

395 

6.2 

Fiji 

2 

05/07/87 

03:05:49.1 

46.74E 

139.23E 

430 

6.0 

Sea  of  Okhotsk 

3 

07/08/87 

11:50:14.7 

26.97S 

108.16W 

10 

6.1 

Easter  Is. 

4 

07/14/87 

23:46:03.5 

49.63N 

147.83E 

576 

6.0 

Sea  of  Okhotsk 

5 

07/16/87 

05:46:29.5 

33.06N 

138.10E 

310 

5.3 

Honshu 

6 

08/08/87 

15:48:56.7 

19.02S 

69.99W 

70 

6.4 

Chile 

7 

10/27/87 

21:58:17.0 

28.68S 

62.93W 

605 

6.0 

Argentina 

8 

03/25/88 

19:36:46.4 

62.I5N 

124.18W 

10 

6.1 

NW  Territories 

9 

04/12/88 

23:19:55.5 

17.19S 

72.31W 

33 

6.1 

Peru 

10 

04/25/88 

10:10:33.8 

7.79S 

158.26E 

44 

6.1 

Solomon  Is. 

11 

08/10/88 

04:38:26.5 

10.20S 

160.80E 

38 

6.1 

Solomon  Is. 

12 

08/14/88 

10:56:57.5 

54.61  N 

152.68E 

645 

5.4 

Sea  of  Okhotsk 

13 

09/07/88 

11:53:25.5 

30.32N 

137.37E 

501 

6.0 

Honshu 

14 

10/08/88 

04:46:27.8 

18.69S 

172.52W 

65 

6.6 

Tonga  Is. 
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Table  2 
Data  Coverage 


Event 

Distance 
from  BKS 
(“) 

Backazimuth 
from  BKS 
(“) 

BKS 

MHC 

SAG 

1 

125.5 

271.9 

* 

2 

68.6 

313.3 

4t 

3 

66.0 

166.3 

* 

4 

62.2 

312.9 

* 

5 

77.2 

302.1 

*- 

6 

75.0 

129.3 

♦ 

* 

7 

86.5 

130.9 

* 

8 

24.5 

357.8 

* 

9 

72.2 

129.8 

* 

♦ 

10 

86.6 

257.4 

* 

* 

11 

86.1 

254.0 

* 

% 

12 

57.5 

316.8 

* 

* 

13 

79.3 

300.2 

* 

* 

* 

14 

73.6 

229.4 

♦ 

4c 

4c 
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Table  3 

Inversion  results  for  Q,  and  Q,  at  MHC 

(Qr.Q.) 

0.75  hz 


Event  1.5  hz 


0.375  hz 


Table  4 


Inversion  results  for  Q,  and  Q,  at  BKS 
(0,.0,) 


Event 

1.5  hz 

0.75  hz 

0.375  hz 

6 

- 

129,170 

- 

9 

- 

139,117 

119,54 

10 

463,259 

146,128 

127,66 

11 

289,509 

133,233 

104,42 

13 

142,754 

117.201 

- 

14 

276,270 

16.3.105 

91.63 

Table  5 


Inversion  results  for  Q,  and  Q,  at  SAO 


Event 

1.5  hz 

0.75  hz 

0.375  hz 

9 

- 

150.133 

119,54 

11 

351,267 

113,154 

234,66 

12 

209,106 

143,232 

75,103 

13 

160,236 

106,144 

- 

14 

329,311 

148.157 

98,1 14 
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Appendix 

The  convolution  of  s(t)  with  r(t)  (equations  3  ard  4)  is  straightforward  except  for  the  term 


(Al) 


This  convolution  is  approached  by  using  the  discrete  form  of  the  time  series,  equations  (5)  and  (6), 
for  both  the  random  component  and  the  decaying  exponential.  Let 


h(t)  =  C,la^e^‘‘%t-JAa) 

J’O 


(A2) 


^  — ly  kAl 

g(t)  =  c,  i  b^e  '  6(t-kAi) 

Jk>0 


(A3) 


The  Fourier  transforms  are 


h(a})  =  C,  I 
/=0 


(A4) 


«(a))  =  C,  IV'  ^ 

4=0 


(AS) 


Performing  the  convolution  h(t)*g(t)  by  spectral  multiplication 


hi(0)Hw) 

L\Cr 


.  f  I  T  .uvu 

f  xlAt  -T,-V  -T,2A/ 

+  \a^2^ 


=  L  1'  a.e  b.^.e  c 


n  * 0/  0* 


(A6) 


-  jO- 


Since  the  aj  and  b„.j  are  random  variates  between  1  and  -1  their  product,  c„j,  is  also  a  random  variate 
between  1  and-1  although  the  product  distribution  function  is  not  Gaussian.  A  further  complication 
is  introduced  by  the  decaying  exponentials  in  equations  (A3)  and  (A4)  since  the  statistical  nature 
of  the  Gaussian  noise  series  does  not  remain  stationary  with  time.  Particular  values  of  the  new 
non-stationary  random  time  series  will  depend  on  particular  numerical  realizations  of  aj  and  b^. 
However,  a  heuristic  approach  can  be  taken  to  gain  some  insight  on  the  nature  of  the  envelope 
function  which  is  of  primary  interest  to  this  discussion.  The  envelope  is  generally  defined  by  those 
values  which  are  maximum  or  minimum.  Clearly,  maxima  and  minima  of  the  convolution  cKcur 
when 

=  1  or  -1  for  7=0  to  n  (A7) 

By  inspection  of  equation  (A6),  we  see  that  at  these  time  points 

/l((l))g(0)  _  ^  ^  '  “YJ-V  -y^in  - f\\i  -KitfiA/  .  4  cj  , 

— 7:7; —  -  L  e  e  '  <•  (Ac's ) 

L,Lr 


which  is  just  the  discrete  time  series  version  of 


(A9) 


In  the  ab.sence  of  the  decaying  exponentials,  equation  (A6)can  be  transformed  into  the  time  domain 
by  using  the  .shift  theorem  to  give 

h{t)*g{t)  =  C,CMt)RU)  (AW) 


where  R(t)  is  the  linear  ramp  function  and  we  define  a  new  random  function 

u(t)=  Z  db{t-jAt)  (All) 

;  =  i 

For  small  j  (short  times)  dj  is  a  statistical  variate  composed  of  the  product  of  two  Gaussian  distributed 
variates  with  zero  mean  and  chosen  between  the  range  -1  and  1.  However,  as  j  tends  to  large  N, 
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the  dj  increasingly  become  Gaussian  distributed  through  the  central  limit  theorem  of  the  calculus 
of  probabilities.  The  ramp  function  appears  since  the  variance  of  the  sum  of  the  Gaussian  variates 
contained  in  (A6)  increases  linearly  with  time. 

Based  on  these  arguments,  an  approKimation  to  the  inverse  transform  of  ( A6)  can  be  given  be  finding 
the  convolution  of  envolope  functions  (the  decaying  exponentials)  as  the  new  envelope  function  of 
ult).  Therefore 

-ir,'  -if,» 

h{t)*g(t)  =  C,C,u(t)  if 

Jr  Is 

=  C,CMnte-'<'  if  Y  =  Y,  =  Y,  iA12) 

Several  numerical  experiments  were  performed  to  investigate  the  behavior  of  the  envelope  of  a 
number  of  realizations  of  (Al).  It  was  found  that  (A12)  was  an  adequate  approximation  for  the 
purposes  of  this  paper  for  the  relatively  slow  coda  decays  observed  in  the  data.  Since  C,  and  C,  are 
always  less  than  unity,  the  importance  of  this  term  is  usually  dominated  by  the  other  terms  in  the 
total  convolution,  equations  (7)  and  (8). 
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BKS 


^  Three-Component 
Seismic  Station 


JMHC 


Central 

California 


Figure  1:  Sketch  map  of  coastal  Central  California.  Station  locations  for 

Berkeley  (BKS),  Mount  Hamilton  (MHC)  and  San  Andreas  Observatory 
(SAO)  are  shown.  Major  faults  in  the  area  are  also  shown. 
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BKS  Station,  Tonga  Isionds  Event  10/08/88 


fimc  (sec) 


Figure  2:  An  example  of  three  component  data  recorded  at  BKS  for  the  Tonga 
Islands  event  of  10A)8/88.  Z,  R  and  T  are  the  vertical,  radial  and 
tangential  motions  for  the  P  wave  arrival,  respectively.  Note  the 
large  arrival  on  the  radial  component  annotated  by  the  arrow.  This 
phase  is  azimuthally  dependent  and  is  largest  for  events  of  this 
azimuth.  Note  also  the  growth  of  the  tangential  component  with  time, 
suggesting  the  progressive  scattering  of  the  incident  wave  with  local 
velocity  heterogeneity. 
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SAO  Station,  Sea  of  Okhotsk  Event  08/14/88 


50_  '  '  '  1  '  '  '  '  '  '  1  ' 
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Figure  3:  An  example  of  three  component  data  recorded  at  SAO  for  the  Sea  of 
Okhotsk  event  of  8/14/88.  Same  scheme  as  Figure  2.  Note  that  the 
tangential  component  is  as  large  as  the  vertical  component  of  motion 
immediately  at  the  P  wave  arrival  time.  This  effect  is  frequency 
dependent  and  largely  disappears  at  lower  frequencies.  It  is  also 
azimuthally  dependent  and  is  related  to  heterogeneous  fault  zone 
structure  in  the  upper  5  km  of  the  crust. 
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Qr  for  SAO  Q’’  ^KS 


Qr  for  MHC  Qs  for  MHC 


Figure  4:  Correlograms  or  scatter  plots  showing  the  correlation  of  Qr  and  Qs 
measurements  from  the  codas  of  events  which  were  simultaneously 
recorded  at  BKS  and  MHC  (top)  and  SAO  and  MHC  (bottom).  The 
different  symbols  denote  Q  parameters  determined  for  three  different 
bandpass  filtered  versions  of  the  observed  coda. 
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Figure  5:  Vertical  component  of  the  P-wave  for  the  9/07/88  Honshu  event 
recorded  at  each  station.  Absolute  amplitudes  in  digital  counts  are 
displayed.  The  baseline  of  each  trace  has  been  ifted  for  dispb  y 
purposes.  Note  that  the  BKS  P-wave  is  about  t\  the  ampli:  :dc  of 
the  P-wave  at  MHC  and  three  times  the  amplituu  of  the  P-wave  at  SAO 
due  to  focussing  effects  under  the  receivers.  Coda  amplitude  levels, 
however,  are  about  the  same. 
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Abstract 


Numerical  simulations  of  teleseismic  wave  propagation  in  a  heterogeneous  layer 
over  a  homogeneous  half-space  are  conducted  to  further  our  understanding  of 
teleseismic  coda  generation  with  the  hope  of  shedding  light  on  the  generation 
of  complex  teleseismic  particle  motions.  Acoustic  and  elastic  finite  difference 
synthetics  are  generated  for  more  than  150  different  layer  over  a  half  space 
models.  I'lie  models  vary  in  scattering  layer  thickness  (h),  heterogeneity  cor¬ 
relation  distance(s)  (a,  or  ar,aj),  and  heterogeneity  standard  deviation  (tr). 
The  synthetic  data  are  examined  using  relative  code  intensity  envelopes,  to 
determine  relationships  that  exist  between  the  model  parameters  and  the  level 
and  rate  of  coda  decay,  and  frequency-wavenumber  analysis  to  determine  the 
wave  typos  present  in  the  scattered  wavefield. 

The  level  of  scattered  energy  is  found  to  be  controlled  by  the  ka  value  and 
the  magnitude  of  the  standard  deviation  (cr)  of  the  perturbations.  Coda  levels 
increase  with  ka  from  A’a<l  to  and  decrease  with  increasing  ka  for 

ka>l.  Scattered  energy  levels  are  found  to  always  increase  with  increasing 
<7.  A  more  com|)lex  relationship  is  found  to  exist  between  the  coda  level  and 
decay  rate,  and  variations  in  the  scattering  layer  thickness.  For  models  that 
vary  in  scattering  layer  thickness  alone,  coda  level  and  decay  rate  changes  did 
not  occur  independently  of  the  values  of  ka  and/or  a. 

The  rate  of  coda  decay  is  found  to  be  controlled  by  the  aspect  ratio  (0^-/0^) 
of  the  heterogeneities  (spatial  anisotropy).  The  two  extreme  models  are  those 
with  homogencx)us  plane-layers,  which  exhibit  the  most  rapid  coda  decay  rates, 
and  those  with  spatially  isotropic  heterogeneities  (aj./aj=l),  which  exhibit  the 
slowest  coda  decay  rates. 

Acoustic  and  elastic  models  exhibit  similar  coda  intensity  envelope  char¬ 
acteristics.  Api)arcnt  scattering  attenuation  of  the  direct  pulse  is  found  to  be 
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strongest  for  models  witlj  spatially  isotropic  heterogeneities. 

Frequency-wavenumber  analysis  of  the  elastic  synthetics  reveals  the  mech¬ 
anism  responsible  for  the  generation  of  coda  in  the  synthetic  data.  For  mod¬ 
els  with  spatially  isotropic  heterogeneities,  the  coda  is  composed  largely  of 
low  apparent  velocity  [LAV]  waves  (scattered  S  and  possibly  body-to-surface 
scattered  waves)  while  the  coda  for  models  with  homogeneous  plane-layers  is 
composed  solely  of  vertically  propagating  layer  reverberations.  Decreasing  the 
heterogeneity  aspect  ratio  increases  the  total  amount  of  scattered  LAV  energy 
produced.  The  onset  time  of  the  LAV  energy  is  also  controlled  by  the  het¬ 
erogeneity  aspect  ratio.  For  models  with  spatially  isotropic  heterogeneities, 
the  LAV  energy  appears  immediately  after  the  initial  arrival.  For  models  with 
homogeneous,  plane- layers,  LAV  energy  never  appears. 

(key  words:  receiver  function,  teleseismic  coda,  P-to-S  scattering,  Gaus¬ 
sian  correlation,  spatial  anisotropy,  relative  intensities,  f-k  analysis.) 

(page  heading:  scattering  effects  for  teleseismic  wave  propagation) 
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Introduction 


The  use  of  tclcseismic  waveform  modeling  to  determine  the  location  of  nra- 
jor  discontinuities  in  the  earth’s  lithosphere  was  introduced  over  a  decade 
ago  by  Burdick  and  Langston  (1977).  A  major  advantage  of  this  method, 
often  referred  to  as  receiver  function  modeling,  is  that  teleseismic  data 
from  a  single,  isolated  receiver  can  be  used  to  determine  structural  features 
in  the  area  of  the  station.  For  Burdick  and  Langston’s  procedure  to  be 
valid  we  must  assume  that  a  homogeneous,  plane  layered  structure  exists 
beneath  the  station  (Burdick  and  Langston.  1977;  Owens,  1984).  Unfortu¬ 
nately,  observed  data  cannot  always  be  modeled  using  simple  plane-layered 
models  (Langston,  1977b;1979);  this  problem  becomes  more  severe  when  >is- 
ing  broad  band  or  high-freciucncy  data  which  are  more  sensitive  to  small  scale 
heterogeneities  (Owens,  1984;  Owens  cl  ai,  1984,1987).  Particularly  trou¬ 
blesome  in  the  teleseismic  waveform  modeling  problem  are  complex  F-wave 
particle  motions  which  sometimes  contain  anomalously  large  transverse  com¬ 
ponents  (Langston,  1979;19S9).  Slightly  more  complex  modeling  procedures 
that  utilize  dipping  layers  (Langston,  f977b,1979),  or  curvexJ  interfaces  (Lee, 
1983;  Lee  and  Langston,  1983)  have  been  used  to  explain  some  data,  but  a 
great  deal  of  the  observed  complex  particle  motion  cannot  be  modeled  using 
simple  lechni<iues  (Langston,  1979;1989;  Owens  ct  al.,  1987). 

Attempts  to  help  better  explain  and  understand  the  wave  propagation 
phenomena  respo.Tsible  for  the  observed  lele,seismic  particle  motions  and  tele¬ 
seismic  (oda  began  in  tin'  early  1980*.  I'eleseismic  wave  pro|)agation  and 
scattering  for  complex,  plane-layered  structures  were  investigated  by  Richards 
and  Menke  (1983).  It  was  approximately  at  this  time  that  the  large  body 
of  work  u  1  coda  generation  for  loccil  eartlupiakes  began  to  see  widespread 
application  to  the  teleseismic  problem.  For  the  local  earlhcpiake  ca.se,  the  in- 


terest  in  seismic  coda,  the  scattered  waves  arriving  after  the  2n'rival  of  the 
major  seismic  phases,  began  in  the  late  19C0*  and  early  1970*  with  the  in¬ 
troduction  of  the  single  backscattering  model  (Aki,  1969;  Aki  and  Chouet, 
1975).  Since  that  time,  many  studies,  both  theoretical  and  numerical,  have 
been  conducted  to  help  determine  the  relationship  betwt'en  the  degree  of  het¬ 
erogeneity  of  the  earth’s  crust  and  upper  mantle,  and  the  behavior  of  seis¬ 
mic  coda.  Theoretical  models  include  the  single  backscattering  model  (Aki, 
1969,1980,1982;  Aki  and  Chouet,  1975;  Sato,  1977),  multiple  scattering  mod¬ 
els  (Gao  etal.,  ,  19S3a, 1983b;  Gcio  and  Li,  1990;  Zeng  et  ai,  1990),  diffusion 
models  (Dainty  <:t  ai,  1974;  Nakamura,  1976,1977;  Dainty  and  Toksoz,  1977; 
Kopnichev,  1977;  Toksoz  et  al.,  1988),  energy-flux  models  (Frankel  and  Wen- 
nerberg,  1987),  energy-flux  with  radiative  diffusion  models  (Langston,  1989; 
Korn,  1990),  and  transport  theory  models  (VVu,  1985,1988).  .Numerical  inves¬ 
tigations  into  the  mechanisms  responsible  for  coda  generation  include  those  for 
complex  plane-layered  structures  (Richards  and  Menke,  1983;  Langston,  1989). 
irr«‘gular  surface  and  sul)surface  layers  (Hill  and  Levander,  1984;  Levander  and 
Hill,  1985),  and  2-D  inhomogenc-ous  half  and  whole-spaces  (Frankel  and  Clay¬ 
ton,  1984,1986;  McLaughlin  etal.,  1985;  Dougherty  and  Stephen,  1988;  Korn. 
1990). 

Despite  the  numerous  studies  investigating  seismic  coda,  many  questions 
about  its  origins  and  relation  to  the  heterogencitic's  of  the  earth’s  crust  and 
upper  mantle  remain  unanswered,  particularly  for  the  teleseismic  case.  To 
help  belter  understand  how  various  heterogeneous  structures  affect  teleseis¬ 
mic  coda  and  waveforms  in  general,  we  have  carried  out  a  detailed  numerical 
investigatitni  into  scattering  effects  for  teleseismic  wave  propagation  in  a  2-1). 
ramlondy  heterogein  oils  scattering  l.iyer  u\'er  a  homogeneous  half-space.  More 
than  LIO  models  are  tested  to  determine  how  the  scattering  layer  thickness  (L). 
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the  heterogeneity  correlation  distancc(s)  (a,  or  and  the  the  magnitude 

of  the  heterogeneity  standard  deviation  (<t)  affect  teleseismic  coda  generation 
and  apparent  scattering  attenuation. 

Scattering  effects  for  both  acoustic  and  clastic  media  are  studied  by  com¬ 
puting  synthetic  seismograms  using  finite  [FD]  difference  techniques  (Alford 
et  al.,  1974;  Kelly  el  uL,  1976;  Fraiikel  and  Clayton,  1984;  Frankel  and  Clay¬ 
ton,  19SG).  To  simulate  teleseismic  wave  propagation,  waveforms  are  computed 
using  a  vertically  propagating,  broad-band  (1-second  Gaussian),  plane  com- 
pressionai  wave  (P-wave)  as  the  source.  This  incident  pulse  reflects  from  the 
frec-surface  at  the  top  of  the  model,  and  is  absorbed  at  the  bottom  to  sim¬ 
ulate  its  passage  back  into  the  mantle.  Receivers  locali'd  every  1  km  record 
the  data  at  the  frcc-surfacc.  The  main  questions  for  this  study  are:  (1)  Can 
the  amplitude  and  rate  of  decay  of  the  seismic  coda  be  used  to  determine  the 
degree  of  heterogeneity  of  the  models?  (2)  What  are  the  sources  of  the  scat¬ 
tered  energy?  and  (.3)  are  the  acoustic  approximations  and  currently  popular 
theoretical  models  applicable  to  the  clastic  teleseismic  problem? 

The  synthetic  waveforms  are  analyzed  using  two  bcisic  methods;  1)  a  com¬ 
parison  of  relative  coda  intensity  envelopes  (Aki  and  Chouet,  1975;  Frankel  and 
Wennerberg,  1987;  Langston,  1989),  and  2)  frequency- wavenumber  [f-k]  anal¬ 
ysis.  Relative  coda  intensity  envelopes  are  computed  for  several  pa.ss-bands 
for  each  model,  and  compared  in  an  attempt  to  reveal  any  relationships  that 
exist  between  the  various  model  parameters  and  the  coda  level  and/or  rate  of 
decay.  Frequency-wavenumber  analysis  is  conducted  in  an  attempt  to  identify 
the  scattering  sources  and  wave  types  present  in  the  scattered  wavefield. 


-65- 


Numerical  Modeling 


Models  that  simulate  a  randomly  heterogeneous  layer  over  a  homogeneous  half 
space  were  examined  in  this  study.  Finite  dilTercnce  numerical  techniques  were 
used  because  they  [numerical  techniques]  are,  at  present,  the  only  techniques 
available  to  compute  waveforms  for  wave  propagation  through  a  media  that 
varies  randomly  both  laterally  and  vertically. 

A  complete  statistical  description  of  a  randomly  varying  media  can  be 
obtained  by  specifying  (Ishimaru,  1978);  1)  a  correlation  function,  2)  a  corre¬ 
lation  distance  (a)and,  3)  a  standard  deviation  of  the  perturbations  (it).  To 
accurately  describe  wave  propagation  through  a  randomly  varying  medium, 
a  wave’s  propagation  distance  through  the  medium  must  also  be  know.  This 
propagation  distance  is  often  given  in  terms  of  the  dimensionless  normalized 
propagation  distance,  Lfa  (Chernov,  1960;  Ishimaru,  1978). 

With  lliest;  model  variables  in  mind,  we  generated  a  combined  total  of  more 
than  150  acoustic  and  elastic  models  that  varied  in  correlation  dist.ance(s)  (a, 
or  ax,aj),  percent  standard  deviation  (a),  and  scattering  layer  thickness  (L). 
The  three  basic  types  of  models  used  are:  1)  1-D  models  with  randomly  vary¬ 
ing,  homogeneous  plane-layers,  2)  2-0  models  with  randomly  distributed,  spa¬ 
tially  isotropic  lieterogeiK’ities,  and  3)  2-D  models  with  randomly  distributed, 
spatially  anisotropic  heterogeneities.  Model  parameters  were  varied  in  a  sys¬ 
tematic  manner  in  order  to  span  a  large  range  of  normalized  wavenumbers, 
normalized  projjagation  distances,  and  standard  deviations.  Perturbations 
were  in  velocity  only;  the  density  was  2.5  g/cm^  for  both  the  layer  and  the 
half-space  for  all  inoch-ls.  The  acoustic  models  had  a  mean  layer  velocity  of 
5  km/sec.  I'lie  elastic  models  were  for  Poisson  solids  (Poisson’s  ratio=0.25) 
with  mean  layer  velocities  of  6  and  3. 16  km/ser  for  the  P-  and  .‘s- waves,  respec¬ 
tively.  All  models  have  half-space  acoustic//’  velc;cities  of  8  km/sec  (5-wave 


half-space  velocity =4. C2  kin/sec). 

Each  s(;t  of  tiiodels  was  gcncraU'd  using  every  possible  periiuitalion  of  the 
following  parameters: 

•  scattering  layer  thickness  —  L  was  15,  30  or  45  km. 

•  standard  deviation  of  the  velocity  perturbations  —  cr  was  5,  10,  15  or 
sometimes  20%  of  the  mean  layer  velocity. 

•  correlation  distance(s) 

o  models  with  spatially  isotropic  heterogeneities  —  a  was  1,  2.5,  5, 
7.5  or  10  km. 

o  models  with  spatially  anisotropic  heterogeneities  —  was  1,  2.5 
or  5  km  for  an  of  10,  25,  40,  50,  80  or  125  km. 


Spatial  velocity  variations  were  computed  using  a  Gaussian  correlation 
function  (Frankel  and  Clayton,  1986) 


(‘‘r  +  ^-(klal  +  klal)/A 

2 


spatial  domain(x,z)  wavenumber  domain(kx, k*). 


The  decision  to  use  a  Gaussian  correlation  function  was  motivated  by  the 
availability  of  theoretical  and  numerical  results  from  previous  studies. 

Figures  1  and  figure  2  show  examples  of  models  with  spatially  isotropic  and 
anisotropic  heterogeneities,  respectively  (please  note  the  dilTc’renc*'  in  the  scale 
of  the  two  figures).  These  figures  show  the  magnitude  of  the  perturbations  only. 
These  perturbations  are  scaled  and  superimposed  on  the  mean  layer  velocities 
to  produce  the  actual  models  used  in  the  finite  difference  simulations. 
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Finite  Difference  Modeling 


Finite  difference  modeling  techniques,  accuracy  and  stability  checks,  and  bound¬ 
ary  conditions  have  been  discussed  in  many  papers  (Alterman  and  Karal,  1968; 
Boore,  1972a,b;  Landers  and  Claerbout,  1972;  Alford  et  al.,  1974;  Alterman 
and  Loewenthal,  1974;  Smith,  1974;  Ilan  et  al.,  1975;  Ilan  and  Loewenthal, 
1976;  Kelly  etai,  1976;  Clayton  and  Engquist,  1977;  Reynolds,  1978;  Levan- 
der,  1988;  Stephen;  1988;  Daudt  etal.,  1989).  We  used  fourth-order  accurate, 
inhomogeneous,  displacement  finite  difference  approximations  to  generate  our 
synthetics  waveforms.  Based  on  our  own  accuracy  and  stability  tests,  we 
choose  to  use  one-sided,  stress-free  boundary  conditions  at  the  surface  (Ilan 
and  Loewenthal,  1976),  absorbing  boundary  conditions  at  the  bottom  [of  the 
model]  (Clayton  and  Engquist,  1977),  and  transparent  boundary  conditions  at 
the  sides  (Wagner,  1990).  The  grid-spacing  and  time  increments  were  0.25  km 
and  0.02  secs,  ix'spectively.  These  increments  allowed  frequency  resolution 
to  4  Hz.  No  intrinsic  attenuation  was  included  in  any  of  the  modeling  (i.e. 
Q  =  oo).  The  source  pulse  used  for  all  modeling  was  a  Gaussian  shaped  (broad¬ 
band),  1  second,  vertically  propagating  compressional  plane  wave.  Acoustic 
models  were  512  km  wide,  while  the  [spatially]  isotropic  and  anisotropic  [het¬ 
erogeneity]  elastic  models  were  256  and  430  kin  wide,  respectively.  Each  model 
had  a  15  km  thii  k  homogeneous  half-space,  which  consisted  of  a  homogeneous 
layer  with  an  absorbing  boundary  (Clayton  and  Engquist,  1977).  The  com- 
puU-r  lime  reriuired  to  lonijilete  the  models  used  in  this  study  (using  indi- 
viilually,  but  .simultaneously,  a  SUN  3/280,  4/280,  and  TAAC-1  application 
accelerator  board)  can  be  measured  in  months.  The  elastic  models  took,  at  a 
minimum,  an  hour  to  compute  1  sec  of  surface  recorded  data. 


Part  I:  Intensity  Envelope  Analysis 


The  main  goal  of  the  relative  coda  intensity  envelope  analysis  is  to  look  for 
relationships  between  the  size,  shape,  extent  and  magnitude  of  the  hetero¬ 
geneities,  and  the  level  and  rate  of  coda  decay.  The  energy  in  the  scattered 
wavefield  is  compared  to  that  in  the  incident  pulse  by  generating  averaged  in¬ 
tensity  envelopes  computed  by  stacking  tlic  relative  intensity  envelopes  from  10 
receivers  located  at  10  km  intervals.  We  used  a  method  introduced  by  Frankel 
and  Wennerberg  (1987)  to  compute  these  power  sp<'(  tral  densities  with  respect 
to  time. 

Power  spectral  densities,  P{u},t),  are  calculated  using  the  following  proce¬ 
dure  (Figure  3): 

1)  Filter  raw  data  using  two  passes  of  a  two  pole  Butterworth  filter  —  non- 
causal,  zero  phcise,  3  db  down  at  the  corners  (center  frequency  ±0.26  Hz). 
The  acoustic  and  elastic  data  were  examined  using  0.5,  1  and  2  Hz  pass- 
bands. 


2)  Each  band-passed  trace  is  normalized  to  the  intensity  of  its  initial  arrival. 
The  normalization  constant  is  the  area  of  the  squared  amplitude  of  the 
incident  pulse  (~  intensity) 
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where  A/  is  the  duration  of  the  initial  pulse.  These  intensity  envelopes, 
therefore,  show  the  intensity  of  the  coda  relative  to  the  intensity  of  the 
initial  pulse.  Absolute  differences  between  models  were  not  examined. 


3)  The  envelope  of  each  band-passed,  normalized  trace  is  calculated  using 
its  analytic  time  signal. 
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4)  For  each  model  and  each  pass-band,  10  envelopes,  from  stations  spanning 
100  km  of  the  model,  are  calculated  and  stacked  to  give  the  final  plots 
of  log  intensity  vs  time. 

To  avoid  a  deluge  of  figures,  only  a  few  selected  results  are  shown.  These 
figures  used  to  illustrate  these  results  are,  of  course,  selected  to  clearly  illustrate 
various  points  and  may  exhibit  certain  characteristics  more  clearly  than  results 
from  the  analysis  of  similar  models,  but  attributes  of  a  single  model  will  not  be 
used  as  evidence  for  phenomena  that  do  not  appear  consistently.  Because  the 
acoustic  and  elastic  envelopes  exhibit  similar  characteristics  for  most  cases, 
sometimes  only  one  figure  [acoustic  or  elastic]  is  shown  to  make  a  point. 

Plane-Layered  Models 

Results  obtained  using  1-D  models  will  not  be  discussed  because  they  were 
unable  to  produce  realistic  coda  decay  rates,  particularly  for  the  elastic  case 
where  no  horizontal  [displacement]  coda  can  be  generated  without  introducing 
lateral  media  perturbations  (for  a  vertically  propagating  pulse). 

Models  with  Isotropic  Heterogeneities 

riie  following  Ilirce  si-<  tiuns  discuss  tin’  results  obtaiueil  from  2-1)  models 
with  isotropically  shaped  heterogeneities.  Again,  the  three  model  parameters 
that  vary  are  the  magnitude  of  the  heterogeneity  variance,  the  heterogeneity 
correlation  distance,  and  the  thickness  of  the  heterogeneous  layer.  Most  of 
these  results  are  rather  simple  and  agree  with  results  from  previous  theoretical 
and  nunw'rical  studies. 


Changes  in  Correlation  Distance 

The  heterogeneity  correlation  distance  is  directly  proportional  to  the  size  of 
the  heterogeneities,  and  inversely  proportional  to  their  number  density  per 
unit  area,  do  determine  the  ka  value  it  is  important  to  r<’member  that  the 
mean  layer  velocities  are  5  km/sec,  6  km/sec  and  3.C4  km/sec  for  the  acoustic, 
and  P  and  S  velocities,  respc<  tively.  The  pass-bands  examin.ed  are  0.5,  1  and 
2  Hz. 

A  simple  relationship  exists  between  the  scattcTed  energy  level  and  the 
ka  value.  The  level  of  scattered  energy  increases  with  ka  for  ka<\,  reaches 
a  maximum  near  ~  1  (where  the  wavelength  is  equal  to  the  size  of  the 
heterogeneity),  and  begins  to  decrea.se  with  increasing  ka  for  A'a>l.  This 
simple  relationship  holds  regardless  of  the  layer  thickness  and/or  the  standard 
deviation  of  the  perturbations.  Figure  1  shows  the  relative  coda  intensity 
[envelo[)es)  for  ka<l  and  ku>].  Figure  5  shows  the  vcTtical  and  horizontal 
displacement  envelopes  for  elastic  models  that  varied  in  correlation  distance 
alone;  k[’a<\,  where  kp  is  the  P-wave  wavenumb<'r,  for  all  envelopes  in  this 
figure. 

While  the  effects  of  changing  the  correlation  distance  are  clear,  equally  clear 
and  just  as  im[)ortant  are  coda  intensity  envelope  <  haracteristics  that  are  not 
effected  by  cf  'nges  in  the  correlation  distance;  namely,  the  rate  of  coda  decay. 
Changes  in  correlation  distance  affect  the  level  of  scattered  energy  according 
to  a  simple  ka  relation,  but  have  no  noticeable  affect  on  the  rate  of  coda  decay. 
This  is  ev  ident  in  all  of  the  figures  and  was  true  for  all  of  the  models  we  tc'sted. 

Changes  in  Standard  Deviation 

Changing  the  standard  deviation  of  the  perturbations  had  a  very  simple  effect; 
incnuLsiiig  the  standard  deviation  increases  the  level  of  scattered  energy.  This 
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effect  holds  regardless  of  the  layer  thickness  and/or  ka  value.  Figure  6  shows 
relative  vertical  and  horizontal  displacement  intensity  envelopes  for  elastic 
models  that  varied  in  percent  standard  deviation. 

As  with  the  models  where  only  the  correlation  distance  was  varied,  changes 
in  the  percent  standard  deviation  effect  the  coda  level  but  have  no  discernible 
effect  on  the  rate  of  coda  decay. 

Changes  in  Scattering  Layer  Thickness 

Effects  caused  by  changes  in  the  scattering  layer  thickness  did  not  follow  any 
simple  relationships.  The  effect  of  scattering  layer  thickness  has  not,  (to  the 
authors’  knowledge)  been  addre.^^sed  in  previous  numerical  simulations,  but 
has  been  postulated  to  plav  a  role  in  coda  levels  and  decay  rates  in  energy-flux 
models  (Langston,  1989;  Korn,  1990). 

The  '■cattering  layer  thickness  has  no  simple,  unambiguous  relationships 
to  either  the  coda  level  or  rate  of  decay.  For  models  with  varying  scattering 
Isy^r  thicknesses  he  coda  level  and  rate  of  decay  seem  to  to  be  related  to 
both  the  ka  value  and  the  standard  deviation.  For  some  models  changing  the 
s<  alu‘ring  layer  thickness  produced  no  marki'd  differenee  in  cither  the  relative 
coda  level  or  rate  of  decay  (Figure  7,  part  a),  but  for  other  models,  there  are 
obvious  changes  in  tlu'  botii  the  level  and  rate  of  coda  decay  (Figure  7,  part 
b).  A  majority  of  the  data  generated  using  models  that  varied  in  scattering 
layer  thickness  alone  showt'd  no  consistent,  systematic  differences  in  either  the 
coda  level  or  rate  of  decay. 

In  some  cases  the  scattering  layer  thickness  does  eflect  the  coda  level  and/or 
rate  of  decay,  but  these  effects  do  not  appear  to  occur  independently  of  the 
value  of  ka  and/or  the  standard  deviation.  Because  of  this  apparently  com¬ 
plex  intcrelationship,  scattering  layer  thickness  changes  are  not  an  appealing 


mechanism  for  explaining  coda  characteristics  because  there  is  no  simple,  un¬ 
ambiguous  cause  and  effect  relationship  between  the  scattering  layer  thickness 
and  either  the  [relative]  coda  level  or  rate  of  decay. 

Varying  Pass-bands 

The  examination  of  different  pass-bands  of  the  same  data  set  yielded  coda 
level/A:a  relationships  similar  to  those  observed  for  the  examination  of  data 
generated  using  models  with  different  correlation  distances.  These  ka  rela¬ 
tionships  were  most  evident  for  and  deteriorated  for  larger  ka.  For 

^a<l,  coda  levels  are  directly  proportional  to  ka,  but  for  A:a>l,  there  are 
no  clear  trends  (within  the  ka  limits  of  this  study).  Thc.se  ka  relations  held 
regardless  of  the  scattering  layer  thickness  and/or  the  standard  deviation.  Fig¬ 
ures  8  and  9  show  the  vertical  and  horizontal  intensity  envelopes  for  elastic 
simulations  where  ka<\  and  ka>l,  respectively. 

And  again,  as  is  the  cases  for  changes  in  correlation  distance  and  percent 
standard  deviation,  the  band-pcissed  data  exhibited  no  noticeable  differences 
in  the  rate  of  coda  decay. 

Models  with  Anisotropic  Heterogeneities 

To  this  point  we  have  found  several  ways  change  the  level  of  scattered  energy, 
but  have  found  no  simple,  systematic,  unambiguous  ways  to  control  the  rate 
of  coda  decay.  This  inability  to  easily  and  systematically  control  the  coda 
av,v.ay  rate  [using  models  with  isotropically  shaped  hetcrogeneitie.s]  motivated 
a  modification  of  our  model  parameterizatioi. /generation  procedure.  The  fol¬ 
lowing  section  discusses  the  relative  intensity  env('lope  relationships  revealed 
from  the  analysis  of  data  generated  using  models  with  anisotropically  shaped 
heterogeneities. 


Approximately  45  acoustic  models  with  anisotropically  shaped  heterogeneities 
were  tested.  The  acoustic  results  allowed  us  to  select  four  elastic  models  that 
would  [best]  characterize  heterogeneity  [spatial]  anisotropy  space.  The  most 
interesting  results  were  those  for  models  with  a  fixed  and  varying  Qx,  The 
two  extreme  models  are  homogeneous,  plane-layered  structures  {ux—oo),  and 
models  with  spatially  isotropic  heterogeneities  (aj.  =  ai),  both  of  which  have 
been  discussed  in  previous  sections.  The  more  [spatially]  anisotropic  models 
mimic  plane-layered  structures  with  layers  of  varying  lateral  extent.  This  type 
of  model  is  intuitively  appealing  because  it  agrees  with  geological  models  that 
contain  layers  that  do  not  have  infinite  lateral  extent,  but  are  truncated  by 
faulting,  folding  or  simple  pinch  outs. 

Figures  10  and  11  show,  respectively,  the  acoustic  and  elastic  intensity 
envelopes  for  (liodels  with  spatially  anisotropic  heterogeneities  with  a  30-km- 
thick  heterogeneous  layer  and  a  10%  standard  deviation.  These  envelopes 
clearly  illustrate  how  the  heterogeneity  aspect  ratio  (spatial  anisotropy)  con¬ 
trols  the  rate  of  coda  d<'cay.  The  decay  rates,  like  the  models  themselves,  are 
bound  by  the  spatially  isotopic  heterogeneity  case  (ai.  =  aj),  which  exhibits 
the  slowest  decay  rate,  and  the  homogeneous,  plane-layer  case  (or^oo),  which 
exhibits  the  most  rapid  decay  rate. 

The  results  using  models  with  anisotropic  het<'rogeneiti«'s  provide  one  of 
the  most  significant  results  of  this  study.  A  simple  cau.se  and  effect  relationship 
<‘xists  between  het<’rog<’n<“ity  ani.sotropy  and  the  rate  of  coda  decay.  The  more 
spatially  isotropic  (  he  hetcTogeneities  (lower  aspect  ratio)  the  slower  the  rate  of 
coda  decay.  Anisotropic  heterogeneities  simulate  more  plane-layered  structures 
and  produce  less  scattered  energy.  These  are  intuitively  appealing  explanations 
and  have  analogies  to  very  heterogeneous,  highly  fractvirefl,  tectonically  active 
ar(;as  which  produce  high  levels  of  scattered  energy,  and  more  stable  cratons 


which  are  presumably  more  plane-layered  and  produce  less  scattering. 


Discussion 

The  results  from  the  relative  coda  intensity  analysis  can  be  best  summarized 
in  two  figures.  Figures  12  and  13  are  composite  plots  showing  relative  coda 
levels  and  decay  rates,  respectively,  for  several  different  models.  Figure  12 
shows  that  coda  levels  are  related  to  the  ka  value  and  the  magnitude  of  the 
standard  deviation  of  the  perturbations.  Scattered  energy  levels  increase  with 
ka  for  ka  <  I  until  A:a~l,  and  then  begin  to  decrease  with  increasing  ka 
for  ka>  i.  Increasing  the  standard  deviation  of  the  perturbations  always 
increases  the  scattered  energy  levels.  These  effects  have  been  observed  in 
previous  numerical  studies  (Frankel  and  Clayton,  198'^, 1986;  Dougherty  and 
Stephen,  1988).  It  is  interesting  to  note  the  similarity  between  the  levels  of 
the  vertical  and  horizontal  components  of  the  elastic  synthetics  (Figure  12). 
This  similarity  shows  that  the  scattered  energy  is  partitioned  equally  between 
the  vertical  and  horizontal  components. 

Coda  decay  rates  can  be  changed  using  different  scattering  layer  thicknc.ss 
for  only  some  models,  but  a  less  ambiguous  and  more  physically  appealing 
way  to  control  the  rate  of  coda  decay  is  to  vary  the  heterogeneity  aspect  ratio. 
Figure  13  shows  that  the  main  factor  controlling  decay  rate  is  the  heterogene¬ 
ity  aspect  ratio.  The  end  members  of  model  space  are  models  with  spatially 
isotropic  heterogeneities,  which  exhibit  the  slowest  decay  rates,  and  the  homo¬ 
geneous,  plane-layered  models  (extreme  spatially  anisotropic  heterogeneities) 
which  exhibit  the  most  rapid  decay  rates. 

Although  not  shown  in  any  of  the  figures,  initial  pul.se  amplitude  attenu¬ 
ation  is  greatest  for  models  with  spatially  isotropic  heterogeneities.  And  the 
initial  pulses  for  synthetics  generated  using  models  with  spatially  isotropic 
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heterogeneities  had  significantly  lower  high  frequency  content. 
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Part  II:  Frequency- Wavenumber  Analysis 


Frequency- wavenumber  [f-k]  analysis  was  used  to  help  identify  scattering 
source  and  the  wave  types  present  in  the  scattered  wavefield.  The  graphi¬ 
cal  display  of  f-k  results  show,  directly,  the  propagation  directions  and  ap¬ 
parent  velocities  of  any  coherent  [plan  wave]  eiu;rgy  traversing  an  array  of 
receivers  (Robinson,  1967;  Gold  and  Rader,  1969;  Beauchamp,  1973,1975; 
Claerbout,  1975;  Ingate  el  ai,  1985;  Kvaerna  and  Ringdal,  1986).  We  used 
one-dimensional  [1-D]  f-k  analysis  on  the  elastic  synthetic  data  (vertical  com¬ 
ponent  only)  to  help  identify  the  scattered  waves  crossing  the  array.  The 
results  were  somewhat  surprising  and  represcMit  an  aspect  of  scattering  that 
has  been  neglected  in  deterministic  and  phenomenological  coda  decay  models, 
as  well  as  many  previous  numerical  simulations  that  rely  on  an  acoustic  wave 
propagation  cissumption. 

To  avoid  any  contamination  of  the  f-k  results  with  model  edge  effects,  the 
f-k  analysis  is  conducted  on  data  that  should  be  void  of  edge  effect  arrivals. 
For  models  with  spatially  isotropic  heterogeneities  (256  km  wide),  the  earliest 
arrival  time  for  any  5-wave  edge  effects  at  the  center  receiver  is  ~37  secs. 
Our  11  stations  array  should,  therefore,  be  free  of  5-wave  edge  effects  for 
~35  secs.  For  models  with  spatially  anisotropic  heterogeneities  (430  km  wide), 
the  earliest  arrival  time  for  5-wave  edge  effects  is  ~60  s(x:s;  none  of  the  analysis 
windows  for  tlie  models  with  anisotropic  heterogeneities  go  beyond  45  secs. 
The  reason  for  the  concern  with  only  5-wavc  edge  effects  will  become  obvious 
shortly. 

Spatially  Isotropic  Heterogeneities 

The  scattered  wavefield  produced  by  the  models  with  [spatially]  isotropic  het¬ 
erogeneities  contains  a  great  deal  of  low  apparent  velocity  [L.AV]  energy  that 


is  in  the  form  of  5-wave3  and/or  body-to-surface-wave  scattered  energy.  This 
LAV  scattered  energy  is  present  in  all  of  the  [elastic]  data  generated  using 
models  with  spatially  isotropic  heterogeneities.  The  heterogeneity  correlation 
distance  appears  to  play  the  major  role  in  the  production  of  this  LAV  energy. 

Figure  14  shows  35  secs  of  the  synthetic  waveforms  generated  using  mod¬ 
els  with  spatially  isotropic  heterogeneities  with  correlation  distances  of  1  and 
10  km  —  the  two  extreme  models  used  in  this  study.  The  correlation  distance 
is  the  only  parameter  that  differs  in  these  two  models;  the  standard  deviation 
and  scattering  layer  thickness  for  both  models  are  10%  and  30  km,  respectively. 

Frequency-wavenumber  analysis  was  used  to  examine  several  different  win¬ 
dows:  (1)  a  window  that  included  all  of  the  data  from  0-35  secs,  (2)  a  5  sec 
window  containing  the  data  between  the  initial  arrival  and  the  first  multiple, 
and  (3)  two  5-.sec  windows  containing  data  arriving  after  the  first  multiple  but 
prior  to  any  5-wave  edge  effects.  The  f-k  spectra  for  these  windows  are  shown 
in  Figure  15.  Several  important  features  are  shown  in  this  figure,  some  that 
are  common  to  both  models  (a  =1  km  and  a  =10  km),  and  several  that  differ. 

The  most  noticeable  and  significant  feature  common  to  both  data  sets  is 
t  •  presence  of  L.\V  energy  immediately  after  the  passage  of  the  initial  pulse. 
This  LAV'  eiKTgy  (onlinues  for  t!ic  rest  of  the  record  and  decays  slowly  with 
time. 

There  are  ahso  several  differences  between  the  two  data  sets.  First,  for  the 
0-35  sec  window,  it  is  obvious  that  the  amount  of  LAV  energy  produced  in 
the  model  with  a  =  10  km  is  mm  li  greater  than  that  produced  in  the  model 
with  «=1  km.  'I’liis  fat  t  is  further  demonstrated  by  the  later  5-sec  windows. 
For  all  of  the  5-sec  windows,  the  a  =10  km  synthetics  contain  more  energy 
(see  max  power  on  the  figures)  than  the  a=l  km  synthetics.  Secondly,  the 
dominant  frequencies  and  wavenuml)ers  of  the  scattered  waves  is,  in  general. 
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much  higher  for  the  small  heterogeneity  model.  This  shows  that  the  size  of 
the  scatterer  is  directly  related  to  the  wavenumber  and  frequency  content  of 
the  scattered  waves. 

Spatially  Anisotropic  Heterogeneities 

The  f-k  results  for  data  produced  using  models  with  spatially  anisotropic  het¬ 
erogeneities  agree  witli  those  from  the  spatially  isotropic  heterogeneity  mod¬ 
els  in  some  aspects,  but  differ  in  another  important  aspect.  Models  with 
anisotropic  heterogeneities  are  capable  of  producing  LAV  scattered  energy, 
but  the  onset  time  of  this  energy  is  delayed;  the  more  anisotropic  the  hetero¬ 
geneities,  the  longer  the  delay.  The  two  (spatially]  anisotropic  heterogeneity 
models  discussed  here  differ  in  horizontal  correlation  distance  only;  is  5  km 
for  both  models,  while  is  25  or  125  km.  The  data  generated  using  these  two 
models  are  shown  in  Figure  16.  The  f-k  analysis  of  this  data  further  empha¬ 
sizes  the  importance  of  the  P-to-S  scattering  mechanisms  in  the  tclcseismic 
coda  problem. 

hrequency-wavenumber  analysis  was  performed  on  several  windows,  some 
identical  and  some  from  similar  sections  of  the  data.  Figure  17  shows  the  f  -k 
spectra  for  the  data  generated  using  the  models  with  [spatially]  anisotropic 
heterogeneities.  It  is  important  to  remember  that  the  low  apparent  velocity 
scattered  energy  generated  by  the  spatially  isotropic  heterogeneity  models  be¬ 
gan  to  arrive  immediately  after  the  incident  pulse  (Figures  14).  Figure  17 
(10-15  and  8-13  sec  windows)  shows  the  f-k  spectra  for  the  data  recorded 
between  the  initial  arrival  and  the  first  reverberation  for  the  two  spatially 
anisotropic  heterogeneity  models.  Unlike  the  data  from  the  spatially  isotropic 
heterogeneity  models,  the  scattercxl  energy  arriving  immediately  after  the  ini¬ 
tial  arrival  cotisists  mainly  of  vertically  propagating  waves  (layer  reverbera- 
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tions)  with  only  a  slight  amount  of  LAV  energy  for  the  0^=25  km  model,  and 
none  for  the  01  =  125  km  model.  This  trend  continues  further  into  the  data 
with  still  no  or  very  little  LAV  energy  up  to  25  secs.  It  is  not  until  after  25  secs 
that  the  LAV  energy  energy  begins  to  appear  for  the  0^=25  km  model.  For  an 
identical  window  in  the  0^=125  km  data,  there  are  still  only  vertically  propa¬ 
gating  waves,  and  the  energy  carried  by  these  waves  is  decaying  very  rapidly. 
Windows  further  into  each  data  set  show  that,  for  the  0^=25  km  model,  the 
familiar  LAV  energy  has  appeared,  while  for  the  ai=\25  km  model,  the  now 
almost  negligible  energy  is  still  mostly  in  the  form  of  vertically  propagating 
layer  reverberations  with  only  a  minuscule  amount  of  LAV  energy. 

Discussion 

The  f-k  analysis  of  the  synthetic  data  provides  several  significant  results.  First 
is  the  importance  of  P-lo-S  and  possibly  body-to-surface-wave  scattering  in 
the  elastic  problem,  and  second  is  the  importance  of  heterogeneity  aspect  ratio 
in  determining  the  magnitude  of  coda  generation.  For  a  model  with  spatially 
isotropic  heterogeneities,  LAV  energy,  in  the  form  of  F-to-S  and/or  body-to- 
surface-wave  scattered  energy,  dominates  the  scattered  wavefield  and  arrives 
iiiunediately  after  the  initial  pulse.  As  the  the  aspect  (aJ■/a^)  ratio  of  the 
heterogeneities  increases,  LAV  scattered  energy  is  less  evident  and  may  not 
appear  till  much  later  in  the  data. 

Heterogeneity  asi)ect  ratio  is  an  intuitively  appealing  way  to  account  for 
coda  decay  rate  differences,  in  addition  to  complex  particle  motions  observed 
in  teleseismic  coda.  Models  with  spatially  isotropic  heterogeneities  are  void 
of  any  plane-layered  strm  ture  and  could  be  used  to  model  highly-fractured, 
tectonically  active  regions  which  are  often  associated  with  higher  levels  of  scat¬ 
tered  energy.  Tiie  spatially  isotropic  heterogeneity  models  arc  also  capable 
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of  producing  significant  amounts  of  horizontal  displacement,  a  characteristic 
observed  in  some  teleseismic  data.  Models  with  spatially  anisotropic  hetero¬ 
geneities  could  be  used  to  model  stable,  plane-layered  regions. 

The  presence  of  P-to-S  and  possibly  body-to-surface-wave  scattered  energy 
is  also  important  in  that  these  phenomena  cannot  be  modeled  with  an  acoustic 
approximation.  The  intensity  envelope  analysis  suggested  that  acoustic  simu¬ 
lations  may  be  a  good  approximation  to  the  elastic  problem;  the  f-k  analysis 
ahs  shown  that  this  is  definitely  not  the  ca.se.  The  presence  of  coverted  energy 
also  suggests  that  most  phenomenological  scattering  models  must  be  further 
refined  to  include  phase  conversions. 
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Summary 

Coda  Levels  and  Decay  Rates 

This  study  has  shown  that  coda  levels  are  controlled  by  both  the  ka  value  and 
standard  deviation  of  the  perturbations,  while  coda  decay  rates  are  controlled 
by  the  «ispect  ratio  of  the  heterogeneities.  Scattered  energy  levels  were  found 
to  increase  for  increasing  ka  from  ia<Cl  to  ko'^l  at  which  point  scattering 
levels  began  to  slowly  decrease  with  increasing  ka.  Scattered  energy  levels 
were  also  directly  proportional  to  the  magnitude  of  ti.e  standard  deviation  of 
the  velocity  perturbations.  Both  of  these  phenomenon  occurred  regardless  of 
the  scattering  layer  thickness,  or  the  standard  deviation  of  the  perturbations. 
These  ka  and  a  effects  occurred  for  both  elastic  and  acoustic  media,  and  to 
an  equal  extent  in  the  vertical  and  horizontal  components  of  the  elastic  data. 

Variations  in  scattering  layer  thickness  had  no  simple,  unambiguous  effects 
on  the  coda.  Level  atul  decay  rate  variations  caused  by  changes  in  the  scat¬ 
tering  layiT  thickness  were  akso  related  to  the  value  of  ka  and  't.  While  (he 
exact  effect  of  layer  thickness  is  not  clear,  it  does  have  an  effect  in  some  cases 
and  should  be  investigated  in  more  detail. 

The  inability  to  systematically  control  the  rate  of  coda  decay  motivated  use 
of  models  with  anisotro|)ically  shaj>«?d  heterogeneities.  This  reparameterization 
allowed  us  to  ea.sily  control  coda  decay  rales  and  provides  an  intuitively  appeal¬ 
ing  mechanism  to  explain  coda  generation.  The  [two]  extreme/end  member 
models  are  those  with  spatially  isotropic  heterogeneities,  and  those  with  ho¬ 
mogeneous,  plane-layers.  The  homogeneous  plane-layered  models  exhibited 
the  most  rapid  decay  rates,  while  those  with  spatially  isotrojne  heterogeneities 
«‘xhil)iled  tin’  slowest.  Decay  rales  anywhere  betw<‘en  the  two  extrem<‘s  (an 
be  obtained  by  altering  the  a.spect  ratio  of  the  heterogeneities.  Physically, 
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these  changes  correspond  to  totally  random  structures  with  no  layering  (as¬ 
pect  ratio=I,  spatially  isotropic  heterogeneities),  and  those  with  homogeneous, 
plane-layers  (aspect  ratio=oo;  extreme  spatially  anisotropic  heterogeneities). 
Tectonically,  these  models  can  be  thought  of  as  representations  of  highly  frac¬ 
tured,  tectonically  active  areas,  and  stable,  more  homogeneously  plane  layered 
cratonic  regions. 

Frequency-Wavenumber  Analysis 

The  f-k  analysis  of  the  synthetic  data  showed  that  wide-angle  (low  apparent 
velocity)  P-lo-S  scattering,  and  possibly  body-to-surface-wave  scattering,  play 
a  major  role  in  telescismic  coda  generation.  In  addition  to  the  importance  of 
surface  and/or  P-to-5-wave  scattering  are  effects  a.ssociated  with  the  hetero¬ 
geneity  a.spect  ratio.  Models  with  spatially  isotropic  heterogeneities  generate 
LAV  energy  that  immediately  follows  the  incident  pulse.  The  higher  the  het¬ 
erogeneity  aspect  ratio  {a^/a^)  becomes,  the  longer  is  the  delay  between  the 
passage  of  the  initial  pulse  and  the  onset  of  LAV  scattered  energy.  For  the 
most  extreme  [spatially]  anisotropic  case  (homogeneous  plane-layers),  LAV  en¬ 
ergy  never  appears  and  the  coda  is  composed  solely  of  layer  reverberations.  In 
short,  the  heterogeneity  aspect  ratio  controls  the  rate  of  coda  decay  by  con¬ 
trolling  the  rale  at  which  vertically  propagating  tele.seismic  <'nergy  is  scattered 
into  the  horizontal  direction. 

The  f-k  analysis  also  showed  that  there  is  a  direct  correlation  between 
the  frequency- wavenumber  content  of  the  scattered  waves  and  the  size  of 
the  heterogeneities.  For  models  with  spatially  isotropic  heterogeneities,  the 
larger  heterogeneities  produce  longcr-wavciength  (lower- frequency)  scattered 
energy,  and  ♦he  smaller  heterogeneities  producing  shoi ter-wavelength  (higher- 
frequency)  scattered  energy. 
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Figure  1 :  Spatially  isotropic  velocity  perturbations  generated  using  an 
isotropic  Gaussian  correlation  function  with  correlation  distance 
of  1,  5  and  10  km. 
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Figure  2:  Spatially  anisotropic  velocity  perturbations  generated  using  an 
anisotropic  Gaussian  correlation  function.  All  three  models  use 
3^=5  km  and  a^lO,  25  and  125  km. 
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Relative  Coda  Intensity  Computation 
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Figure  3:  Procedure  used  to  calculate  relative  coda  intensity  envelopes.  The 
envelope  above  the  band-pass  filtered  trace  has  been  offset  slightly 
for  this  figure. 


Changes  in  Correlation  Distance 


Figure  4:  Relative  acoustic  coda  intensity  envelopes  for  models  with  spatially 
isotropic  heterogeneities  that  vary  in  correlation  distance. 

Correlation  distances  and  the  corresponding  line  weights  are  shown 
in  the  upper-right  comer  of  the  figure.  The  pass-band,  scattering 
layer  thickness  and  standard  deviation  are  also  listed  on  the  figure. 
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Changes  in  Correlation  Distance 


time  (secs) 


both  for  kpa<1 

Figure  5:  Relative  elastic  (vertical  and  horizontal  displacement)  coda 
intensity  envelopes  for  models  with  spatially  isotropic 
heterogeneities  that  vary  in  correlation  distance.  Correlation 
distances  and  the  corresponding  line  weights  are  shown  in  the 
upper-right  comer  of  the  figure.  The  pass-band,  scattering 
layer  thickness  and  standard  deviation  and  displacement  component 
are  also  listed  on  the  figure. 


Changes  in  %  Standard  Deviation 
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Figure  6:  Relative  elastic  (vertical  and  horizontal  displacement)  coda 
intensities  for  models  with  spatially  isotropic  heterogeneities  with 
different  percent  standard  deviations.  Percent  standard  deviation  and 
the  corresponding  line  weights  are  shown  in  the  upper-right  comer  of 
the  figure.  The  pass-band,  scattering  layer  thickness,  correlation 
distance  and  displacement  component  are  also  listed  on  the  figure. 


-99- 


Changes  in  Scattering  Layer  Thickness 


time  (secs) 


Figure  7:  Relative  acoustic  coda  intensity  envelopes  for  spatially  isotropic 
heterogeneity  models  that  vary  in  scattering  layer  thickness. 
Scattering  layer  thicknesses  and  the  corresponding  line  weights  are 
shown  in  the  upper-right  comer  of  the  Figure.  The  passband, 
correlation  distance  and  standard  deviation  are  also  listed  on  the 
figure.  Part  (a)  exhibits  no  marked  differences,  while  part  (b) 
exhibits  very  different  decay  rates. 
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Changes  in  Frequency  Pass-Band 
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Figure  8:  Relative  coda  intensity  envelopes  for  different  pass-bands  of 
elastic  data  generated  using  models  with  spatially  isotropic 
heterogeneities.  The  pass-band  and  the  corresponding  line  weights 
are  shown  in  the  upper-right  comer  of  the  figures.  The  correlation 
distance,  standard  deviation,  scattering  layer  thickness  and 
displacement  component  are  also  listed  on  the  Figures. 
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Changes  in  Frequency  Pass-Band 
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Figure  9;  Relative  coda  intensity  envelopes  for  different  pass-bands  of 
elastic  data  generated  using  models  with  spatially  isotropic 
heterogeneities.  The  pass-band  and  the  corresponding  line  weights 
arc  shown  in  the  upper-right  comer  of  the  figures.  The  correlation 
distance,  standard  deviation,  scattering  layer  thickness  and 
displacement  component  are  also  list^  on  the  figures. 
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Changes  in  with  a  Fixed  a^  (5km) 


Figure  10:  Relative  acoustic  coda  intensity  envelopes  for  models  with  spatially 
anisotropic  heterogeneities  varying  in  horizontal  correlation 
distance.  Horizontal  correlation  distances  and  corresponding  line 
weights  are  shown  in  the  upper-right  comer  of  the  figure.  The 
pass-band,  scattering  layer  thickness  and  standard  deviation  are  also 
listed  on  the  figure. 
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Changes  in  with  a  Fixed  a^  (5km) 
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Figure  1 1 :  Relative  elastic  coda  intensity  envelo|3es  for  models  with  spatially 
anisotropic  heterogeneities  varying  in  horizontal  correlation 
distance.  Horizontal  correlation  distances  and  corresponding  line 
weights  are  shown  in  the  upper-right  comer  of  the  figure.  The 
pass-band,  scattering  layer  thickness,  standard  deviation  and 
displacement  component  are  also  listed  on  the  figure. 
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Figure  12:  Relative  coda  levels  for  a  large  range  of  models.  The  spatially 
isotropic  heterogeneity  plots  are  for  1  Hz  data  for  a  30  km  thick 
scattering  layer.  The  anisotropic  data  are  1  Hz  pass-band  data  for 
a  model  with  10%  standard  deviation  and  a  30  km  thick  scattering 
layer;  ka  is  the  normalized  wavenumber,  ax  and  azare  the 
correlation  distances  in  x  and  z,  respectively  (in  km),  and  is 
the  standard  deviation  in  percent  of  the  mean. 
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Figure  13;  Relative  coda  decay  rates  for  a  large  range  of  models.  The 
spatially  isotropic  heterogeneity  plots  are  for  1  Hz  data  for  a  30 
km  thick  scattering  layer.  The  anisotropic  data  are  1  Hz  pass-band 
data  for  a  model  with  10%  standard  deviation  and  a  30  km  thick 
scattering  layer;  ka  is  the  normalized  wavenumber,  axand  a^are  the 
correlation  distances  in  x  and  z,  respectively  (in  km),  and  is 
the  standard  deviation  in  percent  of  the  mean. 
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Figure  14:  Vertical  and  horizontal  displacements  generated  using  models  with 
spatially  isotropic  heterogeneities  with  correlation  distances  of 
1  km  (top)  and  10  km  (bottom). 
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Figure  15:  Frequency- wavenumber  spectra  for  the  data  generated  using  models 
with  spatially  isotropic  heterogeneities  (Figure  14).  Window 
durations,  apparent  velocities  and  estimat^  maximum  power  contents 
are  shown  with  each  spectra. 
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Figure  17:  Frequeiicy-wavenumber  spectra  for  the  data  generated  using  models 
with  spatially  anisotropic  heterogeneities  (Figure  16).  Window 
durations,  apparent  velocities  and  estimated  maximum  power  contents 
are  shown  with  each  spectra. 
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Abstract 

Teleseismic  data  recorded  at  the  NORESS  seismic  array  are  examined  to  de¬ 
termine  the  degree  and  scale  of  lithospheric  heterogeneities  beneath  the  array, 
and  to  locate  local  scattering  sources.  The  data  are  examined  using  relative 
coda  intensity  envelopes,  frequency-wavenumber  spectra,  and  particle  motion 
polarizations.  The  coda  intensity  envelopes  do  not  exhibit  any  clear  ka  trends; 
the  0.5  Hz  pass-band  data  does,  however,  contain  higher  energy  levels  and  a 
high  energy  packet  ~15  seconds  after  the  first  arrival.  Frequency- wavenumber 
analysis  shows  that  most  of  the  teleseismic  coda  is  vertically  propagating.  The 
window  centered  near  15  seconds,  however,  contains  energy  traveling  to  the 
north-east  at  a  low  apparent  velocity.  Particle  motions  in  this  window  follow 
a  retrogra<l<‘  <’lli|)tical  pattern  characteristic  of  Rayleigh  surface  waves.  I'hc 
scattering  source  presumed  to  be  responsible  for  the  observed  body-to-surface 
wave  scattered  energy  is  Lake  Mjosa  located  approximately  50  km  west-south¬ 


west  of  the  array. 


Introduction 


Teleseismic  data  recorded  at  the  NORESS  array  are  examined  to  determine 
the  degree  of  heterogeneity  of  the  lithosphere  beneath  the  array  and  to  lo¬ 
cate  local  scattering  sources.  The  data  examiiuHl  in  tliis  study  are  from  de<‘p 
focus  events  (>350  km)  located  at  angular  distances  between  135°  and  140°. 
These  events  were  selected  for  several  reasons;  (1)  to  ensure  a  nearly  vertical 
angle  of  incidence  at  the  array,  (2)  to  ensure  that  the  analysis  window  is  tis 
free  as  possible  of  major  reflected,  refracted  and  converted  phases  that  would 
[immediately]  follow  the  initial  arrival  from  a  shallow  event,  and  (3)  because 
deep  focus  events  are  thought  to  have  [relatively]  simple  source  functions.  The 
goal  was  to  obtain  teleseismic  coda  that  was  generated  by  a  single,  vertically 
incident  source  pulse.  We  analyzed  the  data  by  comparing  relative  coda  inten¬ 
sity  envelopes  (Frankel  and  Wennerberg,  1987;  Langston,  1989;  Wagner,  1991; 
Wagner  and  Langston,  1991),  frequency- wavenumber  spectra  (Robinson,  1967; 
Gold  and  Rader,  1969;  Capon,  1969;  Beauchamp,  1973,1975;  Claerbout,  1975; 
Ingate  et  ai,  1985;  Nawab  et  al.,  1985;  Kvaerna  and  Ringdal,  1986;  Gupta 
etal.,  1990a, b),  and  particle  motion  polarizations  (Harkrider,  1964). 

Relative  coda  intensity  envelopes  for  various  pass-bands  are  compared  in 
an  attempt  to  reveal  ka  relationships  that  can  be  used  to  determine  the  cor¬ 
relation  distance  of  the  heterogeneities  in  the  lithosphere  beneath  the  array 
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(A:  is  the  wavenumber  and  a  is  the  correlation  distance  used  to  parameterize 
lithospheric  heterogeneities).  For  a  medium  whose  perturbations  can  be  char¬ 
acterized  by  a  Gaussian  correlation  function,  relative  [teleseismic]  coda  levels 
axe  directly  proportional  to  ka  when  A:a<l,  reach  a  maximum  when  and 

are  inversely  proportional  to  ka  when  A:a>l  (Frankel  and  Wennerberg,  1987; 
Langston,  1989;  Wagner,  1991;  Wagner  and  Langston,  1991).  In  addition,  the 
relative  [teleseismic]  coda  decay  rates  and  the  station  to  station  signal  coher¬ 
ence  increase  with  the  aspect  ratio  of  the  heterogeneities  (ax/aj)  (Wagner, 
1991;  Wagner  and  Langston,  1991). 

Array  data  is  used  because  wave  propagation  directions  and  apparent  ve¬ 
locities  can  be  easily  determined  using  f-k  analysis.  In  a  study  similar  to 
this  one,  Gupta  etai  (1990a,b)  used  f-k  analysis  on  data  generated  by  under¬ 
ground  nuclear  explosions  in  both  the  US  and  Soviet  Union  and  identified  a  low 
apparent  velocity  arrival  approaching  form  the  south-west  as  body-to-surface 
wave  .scattered  <'nergy  produced  by  the  topographic  perturbation  associated 
with  Lake  Mjosa  and  tli<‘  adjacent  mountains.  In  addition  to  using  f-k  analysis 
to  delineate  wave  types,  we  have  used  particle  motion  polarizations  to  further 
discriminate  between  body  and  surface  waves. 
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The  NORESS  Array 


The  NORESS  seismic  array  is  a  25  component  array  located  in  southern  Nor¬ 
way  approximately  78  km  NNE  of  Oslo  (60.74°  N,  1 1 .54°  E).  Data  used  in  this 
study  are  from  the  25  vertical,  short-period  stations,  and  the  4  short-period 
3-component  stations.  The  stations  are  arranged  in  concentric  rings  with  a 
maximum  array  diameter  of  3  km.  Figure  1  shows  the  location  of  the  NORESS 
array  and  the  array  geometry. 

Structurally,  the  NORESS  array  is  located  outside  the  Oslo  graben  in  a  re¬ 
gion  where  the  Moho  dips  slightly  to  the  southeast  (Kanestrom  and  Ilaugland, 
1969;  Sellevoll  and  Warrick,  1971).  Of  particular  interest  for  this  study  is  the 
presence  of  Lake  Mjosa  and  the  adjacent  mountains.  Together  they  account 
for  about  1 ,000  meters  of  topographic  relief.  The  lake  is  located  approximately 
50  km  south-west  of  the  array. 

Event  Descriptions 

Data  for  10  events  recorded  at  the  NORESS  array  were  obtained  from  the 
Center  for  Seismic  Studies  [CSS].  All  events  had  source  depths  greater  than 
300  km  and  were  located  at  ranges  beyond  135°.  Deep  focus  events  were 
selected  to  avoid  complexities  added  by  earth  structure  near  the  source,  and 
because  deep  focus  events  are  thought  to  have  [relatively]  simple,  impulsive 
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mechanisms.  Events  at  great  distances  were  selected  to  insure  that  first  arrivals 
had  nearly  vertical  angles  of  incidence  at  the  array.  Deep  focus,  distant  sources 
were  also  selected  to,  as  best  as  possible,  provide  a  fairly  large  window  that 
contained  only  teleseismic  coda.  To  insure  that  the  coda  do  not  contain  other 
major  phases,  the  hypocenter  locations  and  corresponding  travel-time-table- 
predicted  arrivals  were  carefully  examined  for  each  event. 

Of  the  10  events  examined,  only  3  were  analyzed  in  great  detail  Other 
events  were  rejected  based  on  low  signal  to  noise  ratios;  missing  stations; 
complicated,  non-impulsive  initial  arrivals;  or  because  of  the  presence  of  several 
large-amplitude  arrivals  within  the  analysis  window.  All  three  events  occurred 
in  the  Fiji  Island/Tonga  Trench  region  and  are  presumed  to  be  associated  with 
subduction  of  the  Pacific  plate  under  the  Australia-lndia  plate. 

Table  1  shows  the  dates,  source  depths,  source  latitude  and  longitude,  body 
wave  magnitudes,  and  the  first  and  second  travel-time-table-predicted  arrivals 
for  the  thr<'e  <-vents  analyzed  in  detail.  For  l)revity  only  the  re'sults  from  the 
2/10/87  <’veiit  will  be  pr<rsented  in  detail.  Analysis  of  the  other  two  events 
yielded  nearly  identical  results. 

February  10,  1987  Event 

The  observed  arrivals  for  this  event  do  not  coincide  with  travel-time-table  [jre- 
dictions  for  an  event  with  a  range  of  138.3°  and  a  depth  of  39.5  km.  the  USGS 
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Preliminary  Determination  of  Epicenters  [PDE]  location.  In  fact,  fairly  sig¬ 
nificant  hypocenter  relocation  is  required  to  match  the  observed  arrivals  with 
those  predicted  by  travel- time- tables.  If  the  IJSGS  source  depth  is  correct 
(395  km),  the  epicenter  would  have  to  be  moved  more  than  150  km  to  pro¬ 
duce  ^u:rival  time  separations  observed  in  the  data.  The  similarity  of  the  data 
from  this  and  the  4/1/86  event  suggests  that  the  mislocation  may  be  mostly 
in  depth.  The  observed  arrivals  can  be  matched  with  a  depth  change  and  very 
little  or  no  change  in  the  event’s  epicenter.  Regardless  of  the  mislocation  prob¬ 
lem,  no  travel-time-table  predicted  phases  arrive  within  the  analysis  window. 
The  first  arrival  is  proceeded  by  low-amplitude  energy  which  is  only  slightly 
above  the  noise  level;  this  precursory  energy  was  not  included  in  the  analysis 
window.  Figure  2  shows  data  from  the  four  3-component  stations. 

Coda  Envelope  Analysis 

Analysis  of  the  coda  envelopes  in  hopes  of  being  able  to  determine  the  scale 
and  degRv  of  heterogeneity  beneath  the  array  wa.s  not  as  fruitful  as  hoped. 
The  data  were  examined  by  comparing  relative  coda  intensity  envelopes  for 
several  fre<itiency  pass- bands  in  an  attempt  to  reveal  simple  ka  relationships  {k 
is  the  wavenumber  and  a  is  the  correlation  distance  used  to  parameterize  the 
lithospheric  heterogeneities).  A  Rutterworth  filter  was  used  for  0.5.  1,  1.5  and 
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2  Hz  pass-bands.  The  filter  is  3  db  down  at  the  corner  frequencies  which  were 
at  ±0.25  Hz  of  the  frequency  of  interest.  Unfortunately,  ka  relations  revealed 
by  comparing  various  pass-bands  showed  a  limited  range  of  applicability  in 
synthetic  data  (Wagner,  1991;  Wagner  and  Langston,  1991).  Synthetics  for 
teleseismic  wave  propagation  through  a  randomly  heterogeneous  media  (Gaus¬ 
sian  correlation  function)  with  a  correlation  distance  a,  the  level  of  scattered 
energy  W2is  found  to  be  proportional  to  ka  when  /ba<l,  but  trends  were  not 
well  defined  for  ka>]  (for  the  frequencies  examined). 

The  data  used  to  analyze  the  coda  envelopes  are  from  the  four  short-period, 
3-component  stations.  The  vertical  and  two  horizontal  intensity  envelopes  are 
shown  in  Figure  3.  The  coda  level  and  decay  rates  for  the  1 ,  1.5  and  2  Hz  pass- 
bands  are  fairly  similar,  while  the  0.5  Hz  pass-band  has  a  noticeably  higher 
level  and  a  lower  decay  rate.  The  0.5  Hz  pass-band  is  also  anomalous  in  that 
it  contains  a  second  energy  packet  whose  magnitude  is  equal  to  that  of  the 
initial  arrival.  This  second,  large-amplitude  packet  of  energy  suggests  that 
there  may  be  a  separate  scattering  source  responsible  for  .some  of  the  energy 
at  0.5  Hz.  The  horizontal  0.5  Hz  pas.s-band  data  also  contain  this  anomalously 
high  energy  packet  which  is  greatest  on  the  east-west  component. 
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Frequency- Wavenumber  Analysis 


Frequency-wavenumber  analysis  was  performed  on  the  NORESS  data  to  iden¬ 
tify  the  wave  types  crossing  the  array,  and  to  locate  local  scattering  sources. 
Of  particular  interest  was  the  large  energy  packet  in  the  0.5  Hz  pass-band 
data  arriving  approximately  15  seconds  after  the  initial  arrival.  Frequency- 
wavenumber  analysis  was  conducted  on  the  vertical  component  only. 

Several  different  f-k  analyses  methods  are  used:  broad  band  wavenumber, 
narrow-band  wavenumber,  and  frequency- wavenumber  (see  Figures  4,  5,  and 
6).  Each  of  these  methods  provides  slightly  different  information.  The  broad¬ 
band  wavenumber  spectra  show  power  as  a  function  of  kr  {k^  =  kl  +  k^) 
and  6  (plan  view  angle  of  incidence)  for  no  specific  frequency.  The  narrow- 
band  wavenumber  spectra  show  power  at  a  specific  frequency  as  a  function  of 
kx  and  ky.  The  frequency-wavenumber  spectra  show  power  as  a  function  of 
frequency  and  one  wavenumber/spatial  dimension.  For  this  1-D  [frequency- 
wavenumber]  analysis,  the  2-D  NORESS  array  must  be  collapsed  into  a  linear 
array.  The  broad-  and  narrow-band  wavenumber  spectra  are  used  to  determine 
the  azimuth  of  the  line  used  in  this  projection. 

Figure  4  shows  the  frequency-wavenumber  signature  of  the  incident  pulse 
(event  2/10/87)  using  the  three  dilTcrent  f-k  analysis  methods  (broad  band 
wavenunilKT,  narrow  band  wavenumber,  and  frequency-wavenumber).  All  three 
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spectra  show  that  the  incident  pulse  is  traveling  vertically  (i.e.  infinite  appar¬ 
ent  horizontal  phase  velocity).  This  figure  also  illustrates  some  of  the  differ¬ 
ences  in  the  signatures  of  the  three  different  analysis  methods. 

Figure  5  shows  spectra  for  three  windows  taken  further  back  in  the  coda. 
Almost  all  of  the  coda  had  a  nearly  infinite  apparent  velocity  [vertical  angle 
of  incidence]  as  is  shown  in  Figure  5.  This  suggests,  for  the  scale  lengths 
examined  in  this  study,  that  a  plane  layered  or  extremely  high  aspect  ratio 
(ai/a,)  model  is  a  reasonable  approximation  for  earth  structure  deep  into  the 
mantle  beneath  the  array. 

Analysis  of  the  5  second  window  starting  15  seconds  after  the  first  arrival 
revealed  the  presence  of  a  significant  amount  of  energy  traveling  to  the  north¬ 
east  at  a  low  apparent  velocity.  Figure  6  shows  the  three  different  frequency- 
wavenumber  analysis  methods  for  all  three  events  listed  in  Table  1.  Each  of 
these  spectra  clearly  show  a  significant  amount  of  energy  traveling  to  the  north¬ 
east.  The  narrow-band  wavenumber  and  frequency-wavenumber  spectra  show 
that  this  energy  is  traveling  at  approximately  3  km/sec.  I'he  three  narrow- 
band  wavenumber  spectra  in  Figure  6  are  for  different  frequencies  to  show  that 
the  energy  was  not  limited  to  a  narrow  fretpiency  band. 

Figure  7  shows  an  additional  [l-D]  frequency- wavenumber  spectra  for  a 
60  second  window  that  contains  only  coda;  this  figure  shows  that  the  this  low 
piia.se  velocity  energy  carries  a  majority  of  the  coda's  total  energy. 


The  particle  motion  polarization  for  the  data  in  tiie  window  containing 
the  low  phase  velocity  energy  Wtis  examined  to  determine  if  the  energy  was 
propagating  as  a  body  or  surface  wave.  Figure  8  shows  the  particle  motion 
recorded  at  station  C7  for  the  2/10/87  event.  This  figure  shows  the  unpro¬ 
cessed  vertical  component,  and  a  horizontal  component  obtained  by  rotating 
the  observed  north-south  and  east-west  components  to  a  225°  back-azimuth. 
While  the  data  in  this  window  does  not  exhibit  perfect  retrograde  elliptical 
motion  associated  with  a  fundamental  mode  Rayleigh  wave,  it  is  remarkably 
close  considering  no  processing  has  been  done  to  enhance  the  waveform. 

The  measured  apparent  velocity  of  ~3.0  km/sec  and  a  travel  time  of  ap¬ 
proximately  15  seconds  puts  the  scattering  source  responsible  for  the  Rayleigh 
wave  approximately  45  km  away.  Lake  Mjosa  and  the  adjacent  mountains  rep¬ 
resent  approximately  1 ,000  meters  of  topographic  relief  and  are  located  approx¬ 
imately  50  km  to  the  west-south-west  (see  Figure  1 ).  The  structural/topographic 
perturbation  associated  with  the  lake  and  mountains  is  presumed  to  be  the 
source  of  the  body- to-surface- wave  scattered  energy  ob.served  in  all  three  events 
examined  in  this  study.  This  conclusion  is  identical  to  that  reached  in  .several 
previous  studies  (Gupta  ctai,  1990a, b;  Banister  etal.,  1990),  all  of  which  used 
markedly  different  types  of  sources. 


-121- 


Conclusions 


Teleseismic  coda  from  three  deep  focus  earthquakes  were  examined  to  de¬ 
termine  the  degree  and  scale  length  of  the  heterogeneities  in  the  lithosphere 
beneath  the  NORESS  seismic  array,  and  to  search  for  local  scattering  sources. 
Relative  coda  decay  rates  were  similar  for  the  1,  1.5  and  2  Hz  pass-bands,  but 
levels  were  higher  and  decay  rates  lower  for  the  0.5  Hz  pass-t  ..  This  suggests 
that  heterogeneities  in  the  NORESS  region  are  on  the  order  of  12  kilometers 
or  greater.  In  addition  to  higher  levels  and  lower  decay  ratas,  the  0.5  Hz  pass- 
band  data  contains  a  packet  of  energy  that  arrives  approximately  15  seconds 
after  the  first  arrival.  The  energy  contained  in  this  packet  is  equivalent  to  that 
in  the  initial  pulse. 

Frequency-wavenumber  analysis  yielded  two  results.  First,  most  of  the  ob¬ 
served  coda  had  a  nearly  vertical  angle  of  incidence  suggesting  that,  for  the 
scales  examined  in  this  study,  a  plane  layered  structure  is  a  rerisonable  approxi¬ 
mation  deep  into  the  mantle  below  the  array.  Secondly,  the  high  energy  packet 
of  energy  most  readily  observable  on  the  0.5  Hz  pa,ss-band  coda  envelopes 
corresponds  to  a  low  phase  velocity  (~  3  km/sec)  arrival  coming  from  the 
south-west.  All  three  events  examined  for  this  report  produced  this  low  phase 
velocity  energy,  despite  marked  differences  in  their  spectral  content.  Particle 
motion  polarizations  showed  that  this  low  velocity  energy  exhibited  retrograde 
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elliptical  particle  motion  characteristic  of  a  fundamental  mode  Rayleigh  wave. 
The  apparent  velocity  and  arrival  time  of  the  Rayleigh  wave  put  the  scatter¬ 
ing  source  responsible  for  its  generation  about  45  km  to  the  south-west,  the 
approximate  location  of  Lake  Mjosa.  Lake  Mjosa  and  the  adjacent  mountains 
represent  a,  approximately,  1,000  meter  topographic  perturbation  that  is  pre¬ 
sumed  to  be  responsible  for  the  body-to-surface  wave  scattering.  These  results 
concur  with  those  of  Gupta  et  al.  (199a,b)  who  reached  the  same  conclusion 
using  a  markedly  different  type  of  source. 
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Table  1 

Dates  source  depths,  source  latitude  and  longitude,  body  wave  magnitudes,  and 
the  first  and  second  travel-time-table-predicted  arrivals  for  the  three  events 
analyzed  in  this  study 


Event 

z  (km) 

Lat®. 

Long®. 

Mk 

1*'  and  2*^  arrivals 

471/86" 

541 

136.6 

17.9  s 

178.6  w 

5.8 

m 

pP"  @124  secs 

2/10/87 

395-*- 

138.3 

19.5  s 

177.5  w 

6.2 

12a 

pP"  @122  secs 

3/10/88 

617 

139.5 

20.9  s 

178.7  w 

6.0 

■2a 

pP"  @140  secs 
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Figure  1 :  Location  of  the  NORESS  array  and  the  array  geometry.  The  solid 
triangles  represent  vertical  component,  short-period  instrument 
locations.  The  empty  tiangles  denote  the  location  of  three- 
component  sensors. 
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3-Component  Short-Period  NORESS  Data 

,  (2/1 0/87  event) 


1 0  secs 


Figure  2:  Vertical  and  horizontal  short-period  NORESS  data  for  the  2/10/87 
event.  The  station  and  maximum  amplitude  (nanometers)  are  shown 
above  each  trace.  Vertical,  north  and  east  components  are  denoted 
by  z,  n  and  e,  respectively. 


Figure  3:  Stacked  intensity  envelopes  for  various  pass-bands  for  the  vertical 
and  two  horizontal  components  for  the  2/10/87  event.  The  components, 
various  pass-bands  and  corresponding  line  weights  are  shown  in  each 
panel. 
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Figure  4:  Wavenumber  spectra  for  the  0-15  second  window  of  the  vertical 
component  for  the  2/10/87  event. 
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Figure  5:  Wavenumber  spectra  for  randomly  selected  windows  of  the  vertical 
component  for  the  2/10/87  event. 
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Figure  6:  Frequency-wavenumber  spectra  for  the  15-20  second  window  for  three 
events. 
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Figure  7:  One-dimensional  f-k  spectrum  for  the  15-75  second  window  of  the 
vertical  component  for  the  2/10/87  event. 
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Figure  8:  Particle  motion  for  the  unprocessed  vertical,  and  rotated  (to  45 
degrees)  horizontal  components  of  motion  at  station  C7  for  the 
2/10/87  event  The  circle  on  the  particle  motion  diagram  marks  the 
window  start,  and  arrow-heads  point  in  the  direction  of  increasing 
time.  The  exact  analysis  window  is  shown  with  the  data. 
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